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Abstract 

The aim of the work is focused on the generation 
of FEM model by Ansys code, able to effectively 
evaluate the miscuffing errors. For this purpose it 
was measured the circumferential distribution of 
the pressure produced by the cuff on the arm 
during the inflation and deflation phases, 
concerning bladders of longitudinal dimensions 
comprised between 60% and 100% of the 
circumference of the experimental arm. The model 
parameters have been identified by searching 
vessel-occluded conditions at the subject systolic 
blood pressure for the bladder I arm ratio equal to 
80%. The Young's modulus is fixed to 40 kPa and 
the Poisson 's ratio to 0,44. The measured pressure 
distributions for the 60%, 70%, 90% and 100% 
ratios have allowed to verify miscuffing errors 
reported in bibliography by the constructed FEM 
model. 

1. Introduction 

For the last two decades the authors have been 
involved into blood pressure measurements, with 
the development of phenomenological simulators, 
Vallascas (2005) [21], and also with the proposal 
of an innovative application of oscillometric 
method, Vallascas (2012) [22]. This work falls 
within the topics of interest and it was focused on 
the implementation of a finite elements model, 
which is able to highlight miscuffing errors and 
provide capable information in order to use, and 
eventually redesign, the cuff. 
It is well known that the most common mistakes 
during brachial and radial blood pressure 
measurements, come from the adoption of not 
proportionate cuff to the arm circumference 1 . This 
kind of problem is very common into the pediatric 
field, where the patient's arm size has even wider 



1 The phenomenon of miscuffing occurs when the size of 
the bladder, in relation to the arm's circumference is or 
too small, undercuffing, or too large, overcuffing. 



range, of possible values, than on the adult 
population. 

The inner bladder is the most significant 
metrological factor, which must be constrained to 
comply the correct longitudinal and transverse 
dimensions. 

This bibliography field is quite large and 
articulated, mainly including clinical 
measurements, recommendations of regulatory 
Bodies and researches about modeling and 
simulation. Alexander et al. (1977) [1], indicate 
the cuffs selection criteria, until they create a new 
pediatric cuff design. In Reid and Chantler (1982) 
[18], have been reported the corrections made by 
Maxwell, who was one of the first to suggest in 
80% the optimum ratio between the bladder and 
the arm circumference, and to define the 
functional link between measurement errors and 
the quoted ratio in a quantitative way. 
Several authors have highlighted the reading 
difference that it is possible to obtain comparing 
cuffs of different sizes. These results are reported 
in relative terms and don't care about errors in 
absolute terms. Iyriboz er al. (1994) [10], have 
quantified the reading differences obtained using 
12/23 and 15/33 cuffs. Bakx er al. (1997) [3], have 
determined the effects using cuffs having different 
width and height bladder. Marks and Gronch 
(2000) [13], come to the conclusion, based on the 
comparison between the auscultatory and direct 
pressure measurements, that the measurement 
error can be minimized by using a bladder having 
width equal to 46% of the arm circumference. 
Jilek and Stork (2010) [12], had analyzed the 
influence of the bladder's width during the use of 
the radial site. 

The American Heart Association (AHA) (2005) 
[16], sets four cuffs standard sizes, moreover it 
established that "the ideal cuffs should have a 
bladder length of 80%, and a width of at least 40% 
of the arm circumference". A length/width ratio of 
2:1 appear to be optimal, especially regarding the 
insertion problem around the arm that smaller 
ratios might entail. 
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The European Society of Hypertension (ESH) 
(2003) [15], recommend the three cuffs mentioned 
by the British Hypertension Society (BHS) and 
refers to the dimensions standardized by the AHA. 
The numerical modeling works are mostly 
oriented to deepen issues related to the 
oscillometric method. For instance, Foster and 
Turney (1986) [9], Baker et al. (1977) [4], and 
Jeon et al. (2007) [11], based their model on the 
static pressure-volume relation, Ursino e Cristalli 
(1996) [19], propose a whole model with lumped 
parameters, Wang et al. (2001) [23], focus their 
proposal on fuzzy logic, Pinheiro (2008) [17], 
analyzes the signal transduction related issues, 
while Mersich and Jobbagy (2009) [14], had 
developed a set-up in order to determinate the 
cuffs pressure -volume transfer function, also 
analyzing the compression conditions (tight, loose 
and canvas) of themselves on the support, and 
quantifying the measurement errors. 
The authors, Usai and Vallascas (2011) [20], had 
confirmed in a previous memory, by using a 
preliminary and simplified EF mathematical 
model, the data reported by Reid an Chantler 
(2002) [18], for the 12/23 cuff. Furthermore, 
referring to the ratio between transversal and 
longitudinal dimension, they found that the error 
introduced can be considered constant and 
negligible over the 25% of this ratio. The models 
proposed, so far to interpret and simulate 
phenomena of blood pressure measurement, even 
if they have the theoretical formulation, 
sufficiently comprehensive and rigorous in the 
assumptions and interpretation of phenomena, are 
often supported by experimental validation of a 
low level. 

In the present paper the authors have focused their 
attention on the adoption of a method of 
measurement and the development of a 
computational model suitable for the evaluation of 
the errors resulting from the miscuffing and from 
the positioning of the cuff. In this analysis they 
will use pressure distributions on the real human 
arm, previously obtained experimentally. This 
feature requires the availability of a suitable 
measurement set-up and of a 3D "sectorialised" 
model. That permits the application of punctual 
loads on the outer contours of the arm. It appears 
that oriented researches in this direction has not 
ever been proposed in the literature. The modern 
knowledge about the distribution of cuff's 
pressure on the arm are qualitative, and imply a 
uniform distribution equal to the amount of 
pressure inside the cuff. An exception is reported 
by Cristalli, Neuman and Ursino (1994) [7], which 
have conducted pressure measurements on the 



interface between the cuff and the arm, but limited 
to a few points without characterizing the 
distribution on the complete circumference. 
Usually it is ignored the possibility that the cuff 
can be set lobes up onto the cuff, in order to 
locally load the arm, generating non-uniformity of 
pressure. 

2. Setup 

In order to determinate the pressure distribution 
exerted by the cuff on the lateral surface of the 
arm, it was necessary acquiring a suitable set of 
sensors. After it was experimented with a variety 
of available sensors on the market, it was decided 
to use a more reliable set of self-made sensors. 
The sensor is built, Figure 1, with a silicon duct 
with an internal diameter equal to 4 mm, and a 
length equal to 35 ± 1 mm, it is connected to a 
pressure transducer 2 by a smaller diameter rigid 
tube. The transmission fluid is made of distilled 
and degassed water. All the sensors are calibrated 
daily before each measurement cycle. During the 
calibration process are removed the offsets and it 
is calculated the gain, through a comparing 
procedure with an aneroid manometer that derives 
the traceability by a mercury differential U 
manometer. 




Figure 1. Constructive diagram and picture of the 
sensors 

In order to insert the sensors on the measurement 
spot, they have been assembled on a flexible cloth, 
keeping a constant gap one each other; in this way 
they have been able to cover the entire arm 
circumference. 

For not modify the physical distribution of the 
pressures and to contain the sensors, an equivalent 
rigidity material fills the empty space between 
them. During the operational stage the cloth is 
wrapped around the arm, under the cuff, Figure 2. 
The signals have been acquired by a board, using 
NI USB 6251 and Labview 2012 software. 



2 Motorola MPX5050. 
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Figure 2. Operational setup 



The maximum uncertainty of the single 
measurement was estimated to be less than 1 
mmHg. 

3. Measurements 

During the miscuffing analysis it was used a cuff 
12/29, whose bladder, Figure 3, can be extracted 
and one of its ends can be folded back on itself to 
the size of interest, to allow the partializing, 
through the forced insertion of a C aluminium 
profile. 




Figure 3. Bladder used for the partializing of the 
cuff 

Cuff/arm ratios (CC/AC) between 60% and 100%, 
with step equal to 10%, have been analysed. The 
arm circumference, on which the measurement 
were conducted, was measured equal to 260 ± 5 
mm on the cuff insertion site. 
In order to verify the repeatability of the 
phenomenon, the measurements have been 
repeated after the cuff insertion on the arm for 
each configuration. The repeatability, also in 
consideration of the fact that deflation was 
conducted manually, was contained within 5 
mmHg in all the measuring points, Figure 4. The 
picture present on the abscissa the time, expressed 
in seconds and on the ordinate axis the 
corresponding value of the pressure, expressed in 



mmHg. The yellow curve in the center refers to 
the cuff internal pressure during the phases of 
inflation and deflation. 




Figure 4. Example of pressure monitoring 

The Figure 5 shows the diagrams relating to the 
circular distribution of the pressures measured by 
10 sensors for the analyzed configurations. It is 
evident the lack of axial-symmetric distribution in 
all cases. This effect was observed both in 
measurements conducted on human arm than on 
rigid cylinders, with the difference that human 
limb pressure distribution does not undergo 
appreciable variation turning the headset within ± 
20°, while on a rigid support, of course, the 
diagram follows the rotation. From the diagrams 
were derived circumferential pressures, used to 
load the model. 



Figure 5. Diagram of pressure distributions [kPa] 
around the arm 

4. Fern model 

The arm simulation model was developed by 
using Ansys code. It is considered a half of the 
total volume, assuming the model symmetry to the 
YZ plan (cross-section). The software Ansys, 
allows the rebuilding of the sizes distributions of 
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the field, on the whole volume, through post- 
processing operations. 

The arm's tissue has been assumed homogeneous 
and elastic, this because, after the "pressurization" 
accomplished by the cardiac pump, the arm is able 
to resume to the initial configuration. It is also 
wanted to keep the presence of the basilica vein 
inside of the tissue for a better correspondence to a 
real situation. 

4.1 Parameters characterizing the 
materials 

In reference to the cross-sections, quoted in 
Eycleshymer and Schoemaker [8], it has been 
assumed a configuration built with a tissue with a 
circular section, placed in a uniform way around 
the bone and containing the brachial artery and the 
basilica vein, along the same diameter. 
The external bone diameter has been assumed 
equal to 28 mm, the external arm diameter has 
been assumed equal to 108 mm, the artery inner 
diameter has been assumed equal to 4 mm with a 
thickness of 0.6 mm, while the inner diameter of 
the vein has been assumed equal to 7 mm with 
thickness equal to zero. 

In order to further reduce the number of the nodes 
and of the elements, it was considered only the 
bone outer surface, removing the whole volume. 
The rigidity has been imposed to the bone outer 
surface through the CERIG Ansys command. The 
numerical model consists of a total of 38 volumes, 
194 areas, 309 lines, 158 key points, 58284 nodes 
and 39138 solid elements, 92 plus 1 elements, 
Target 170 for the pilot node. The processing time 
was reduced to about 8 minutes. 
The overall arm tissue was modelled as forming a 
set composed of skin, muscle, nerves and fat, all to 
be homogeneous, elastic and isotropic, at this set 
was assigned a modulus of elasticity E = 40 kPa. 
The Young's modulus was chosen, through a 
process of calibration of the model, on the basis of 
experimental values of cuff pressure 
corresponding to the closure of the artery. The 
Young's modulus, is of the same order of 
magnitude than those determined by Zheng on the 
forearm (1999) [24], and that might even be 
overestimated considering the values found for a 
bovine muscle by Whereas Chen (1996) [6]. are 
lower of an order of magnitude. 
The venous vessels' presence has been preserved 
in order to constitute a possible preferential 
breakdown point to facilitate the artery collapse. 
The arterial was considered isotropic and 
characterized by an elastic modulus E= 440 kPa 
and a Poisson's ratio v = 0.45. 



The elastic modulus values are in agreement with 
the results reported in Avolio (1980) [2], and Bank 
(1996) [5]. 

4.2 Constraints and boundary conditions 

As mentioned, the analysed structure concerns 
half of the entire volume of the arm, and it were 
applied to it the following constraints and loads. 
It was imposed a uniform pressure of 100 Pa to the 
discharge zone of the arm, and the pressure 
determined by the experimental distributions to the 
loaded area in correspondence of the cuff. 
It is supposed a pressure of 0.2 kPa within the 
vein, and a diastolic pressure equal to 10.66 kPa 
into the artery. 

The condition of symmetry is applied on the 
transverse surface, which contains the cuff, while 
all the points of the cross section, on the shoulder 
side, were left without constraints. 
It was made the static analysis for small 
deformations. Various tests have been made, and it 
was verified that using the mesh recommended 
level value, equal to 6, the numerical method 
allows to converge to the results in a range 
between 5 and 8 minutes of computation, using a 2 
GHz portable PC with 4 GB of RAM and 220 GB 
of hard drive. 

4.3 Validation 

Since there is an almost unanimous acceptance 
in considering equal to zero the measurement error 
induced by the cuff, when the ratio bladder/arm 
(CC/AC) is around 80%, for the model validation 
is simply verified that with the actual distribution 
of the pressures, determined experimentally, the 
cross section of the artery was occluded, Figure 6. 




o .uiM«) BnS Rnff 



Figure 6. Displacements diagram for CC/AC = 80% 
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For this assessment, given the three-dimensional 
nature of the displacements, it was decided to 
identify the achievement of the occluded vessel 
conditions by means of a visual estimation of the 
diagrams obtained by the numerical program. 

5. Results 

It has been analysed cuff/arm ratios between 
60% and 100% with step 10%. 
Figures 7 and 8 represent the vertical 
displacements diagrams of the structure, relative to 
ratios equal to 60% and 100%, during the working 
stage. 




■0023?4 .WO** .030W ,8»»7 .019745 



Figure 7. Displacements diagrams for CC/AC equal 
to 60% 
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Figure 8. Displacements diagram for CC/AC equal 
to 100% 

Figure 9 represents the summary of the results 
obtained by the FEM model, compared to the 
errors in measuring the systolic pressure, derived 
from the work [16] the results are strictly related 
to the selected 12/29 cuff that is suitable to cover 
100% of the test arm. 



For each the examined situations, the deviation are 
always smaller than 2%, if compared to the 15/30 
cuff. 

The results are, in the opinion of the authors, 
widely acceptable, also considering the uncertainly 
related to the errors that are recruited as a 
reference and which cannot be estimated less than 
a 100 Pa. 
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Figure 9: Measurements errors by the EF model 

The same FEM model has been applied to both the 
determination of the influence of cuff rotation and 
the evaluation of the changes made by forcing the 
cuff on the arm. 

For these experiments were used a commercial 
cuff with bladder 10/24 positioned on the same 
instrumented arm. The outer circumference was of 
290 mm (CC/AC = 82%). 

In the first case are analyzed rotation angles 
comprised between -20° and +20° in steps of 10°. 
The differences fall within ±500 Pa. 




Figure 10. Pressure [kPa] around the arm as a 
function of the cuff length 

In the second analysis, not having a dynamometer 
for measuring the cuffs force on the arm, were 
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taken into consideration three configurations: cuff 
close with as much force as possible, loose cuff of 
10 mm and loose cuff of 20 mm. The related 
circular diagrams are shown in Figure 10. The 
outer curve (length equal to 290 mm) corresponds 
to the maximum forcing, the inner curve (310 mm) 
corresponds to the minimum. The model showed 
no appreciable differences in the situations 
concerned. 

6. Conclusions 

A FEM model, implemented by commercial 
software Ansys and validated by a static analysis 
for small deformations, has been presented. This 
model introduces the possibility of applying 
pressure loads distributed non-uniformly on the 
outer contour of the arm. For the determination of 
the pressure distribution in the interface arm-cuff 
to use as load on the FEM model has been 
prepared an original set up that uses a cuff with a 
bladder-height invariant and variable length 
obtained through an original solution which plans 
to fold back a 'end to the desired length. 
A further test of the model validation was carried 
out to evaluate the error resulting from the use of 
badly dimensioned bladder, in order to quantify 
the phenomenon of miscuffing. 
The results derived from the model developed, 
loaded with the experimental pressure 
distributions, have confirmed, for all the examined 
situations, the diagrams derived from the work of 
Reid & Chantler. 

The FEM model was then applied both to the 
determination of the influence of the rotations of 
the cuff with respect to the location of the artery 
and to the evaluation of the changes made by the 
forcing of the cuff on the arm. 
In the first case were analyzed rotation angles 
comprised between -20° and +20° with step 10°, 
and were found differences in the evaluation of the 
systolic pressure contained within ± 500 Pa. 
In the second analysis, were taken into account 
three different configurations: cuff close with the 
maximum force as possible, cuff loose of 10 mm 
and cuff loose of 20 mm. 

The results showed no appreciable differences for 
all situations considered . The model and the 
method proposed have allowed to highlight and 
confirm knowledge commonly accepted in the 
practice of blood pressure measurement 
concerning a) the optimum dimensional 
relationships, b) the correct positioning of the cuff, 
c) its forcing on the arm and d) the errors 
introduced by the less considered phenomenon of 
miscuffing. In the light of the obtained results the 



authors believe that the level reached by the 
simulation model can afford both to more 
effectively analyze aspects related to the 
measurement of blood pressure and to perform test 
operations of the innovative cuffs that can be 
designed to provide greater flexibility and lower 
measurement errors. 

Therefore the model was validated with 
experimental tests in laboratory. Through 
validation were chosen the most suitable 
parameters to characterize the biological materials 
of the model (blood, bone, artery, vein). 
In the model was performed static structural 
analysis, referring to the condition of equilibrium 
fluid-dynamic that occurs when the artery is 
occluded. 

Therefore if we consider the real nature of multi- 
physics of the problem examined, this work 
constitutes a valid starting point to get also more 
accurate results, applying the FEM multi-physics 
model with a structural analysis that taking 
account of the fluid-dynamic of the real model 
before it occurs artery occlusion. 
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Abstract 

The paper presents a method for mapping the 
temperature distribution of very fast transient 
events (i.e., having a bandwidth of 10 kHz or 
more) by means of a standard infrared camera 
working at 25 Hz frame rate with 320x256 pixels 
full frame. The proposed method is based on 
triggering multiple time-delayed acquisitions 
(MTDA) of the observed thermal phenomenon, 
which has to be reproducible, by means of a very 
precise and stable programmable digital micro- 
controller and by reconstructing the time domain 
IR sequence using the frames acquired at each 
trigger event. The method could find application 
in assessing the reliability of power electronic 
devices and, in particular, to measure 
dynamically the temperature distribution over the 
source metal of Power MOSFETs or IGBTs, 
which feature very fast thermal transients, even in 
the hundreds of microsecond scale and might 
develop local hot spots as a consequence of aging 
or failures. 

1. Introduction 

Measurement of temperature distributions which 
are characterized by very fast transients, below 
the millisecond scale, has always been a difficult 
problem to be addressed and, only in the recent 
past, few experimental techniques have been 
proposed to deal with this challenging task. 

One typical application field, in which recording 
of fast thermal transients is needed, is the 
characterization of power electronics equipment 
[1]. "Operating" temperature has, in fact, 
important consequences for the performance and 
reliability of semiconductor devices. For instance, 
the speed, or maximum operating frequency, of a 
microprocessor typically decreases as the 
temperature increases, and the gain or 
transconductance of a transistor may either 
increase or decrease with increasing temperature, 
depending upon the device type and operating 
conditions. It is also commonly assumed that the 
safety margin or reliability of a semiconductor 



device decreases as the temperature increases [2- 
3]. It is not surprising, then, that a significant 
amount of effort goes into accurately measuring 
the temperature at which devices operate [4-6] . 

Often designers are more interested in achieving 
the complete 2D temperature distribution across a 
well-defined region of the device rather than 
inferring the average temperature of the device 
itself [7]. While in this last case electrical 
methods exploiting the temperature sensitivity of 
specific properties of semiconductor devices can 
be used (e.g., the forward voltage of a pn-junction 
at a constant current is known to vary with 
temperature in a predictable way), detection of 
temperature distribution either requires a single- 
point x-y raster scanning system or the usage of 
an infrared camera. Moreover, whereas electrical 
methods make use of electrical connections which 
are already available for normal device operation 
and are the only type that can be made on fully 
packaged devices [7], measurement of the 
temperature distribution by means of contact or 
non-contact methods always requires a physical 
or optical access to the device surface, which 
therefore must be preliminary opened. 

The thermal excursions that power devices 
undergo in the transient regime can be very large 
and fast. In addition, significant thermal gradients 
can arise on the device surface because of non- 
homogeneous current density distribution, which 
in turn can cause thermal instabilities and device 
failure. Hence, in order to get the map of 
temperature distribution, a measurement method 
with both broad bandwidth and proper spatial 
resolution is needed. The fulfillment of these 
requirements is not trivial. 

Castellazzi et al [1] have presented an 
experimental method based on a polycrystalline 
IR-fibre which conveys the thermal radiation 
emitted by the heated surface within the spectral 
range from 8 to 12 //m onto a photo -detector. The 
system is able to perform single -point transient 
temperature measurements well down to the 
microsecond-range with an estimated spatial 
resolution of about 700 pirn, but both the detector 
and the amplification circuitry needed to be 
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immersed in liquid nitrogen (77 K) in order to 
increase the signal-to-noise ratio of the 
equipment. In addition, a spray coating consisting 
of a diluted carbon-based paint was used to 
uniform the emissivity of the MOSFET surface. 

An improved method was developed by the 
research group of Breglio et al [2, 8-9], who 
present an experimental set-up based on a 
radiometric microscope with a cooled InSb 
single -point sensor. The 2D scanning of the 
measured surface is obtained by means of an x-y 
motorized stage. In this way, by considering the 
microscope spot-size, it was possible to achieve a 
spatial resolution of less than 10 pirn over an 
extended area of 25 x 25 mm 2 , with quite good 
time resolution (less than 10 ms). 

The same research group also developed a 
second system based on thermoreflectance 
mapping, which exploits the thermo -optic effect 
to measure the superficial temperature with a time 
resolution in the order of tenths of nanosecond 
and a temperature resolution of some kelvin, 
along with a much better spatial resolution with 
respect to the previous method [8] . 

An experimental approach similar to that used 
by Breglio et al [2,8] was adopted at the 
University of Messina to develop a reliability 
model for power MOSFETs working in avalanche 
mode [10-11]. 

All the systems described so far rely on single- 
point measurements and hence require a proper 
scanning system in order to reconstruct the 
temperature distribution map over the region of 
interest. Modern infrared cameras with fast focal 
plane array (FPA) would in principle present 
definite advantages compared to single -point x-y 
raster scanning systems. Unfortunately, these 
devices have a common limitation regarding the 
maximum sampling rate at full frame, which is 
usually limited to some tens of hertz. Higher 
frame rates are sometimes possible at the expense 
of a reduced number of pixels. 

Recently, Riccio et al [12] used the equivalent 
time sampling (ETS) concept to measure transient 
temperatures with an equivalent bandwidth of 100 
kHz by means of a standard IR camera, thus 
allowing to overcome the limitations associated 
with single -point measurements as far as the 
observed phenomenon can be repeated 
periodically (as it is the case for power electronic 
devices heated by pulsed currents). 

Here we exploit the same concept to propose an 
improved method that, in addition of allowing full 
frame measurements of replicable transient events 
with high bandwidth (limited only by the 
integration time of the IR camera) and good 
spatial resolution (depending on the optical lens 



used in the experimental set-up), also provides 
automatic compensation of the background 
temperature, thus eliminating systematic errors 
caused by ambient temperature shifts at each 
trigger event. The proposed method makes use of 
a precise and stable programmable digital micro- 
controller used to trigger multiple time -delayed 
acquisitions (MTDA), allowing the time domain 
infrared sequence to be reconstructed based on the 
time-delayed frames acquired. 

To prove the effectiveness of the proposed 
MTDA method, experimental tests have been 
carried out on a planar resistor heated with 
repeated pulsed currents. A practical application 
dealing with a dynamical analysis of the 
temperature distribution over the source metal of 
an IGBT power device is finally presented. 



2. The MTDA method 

A standard cooled IR camera (FLIR SC7000) 
equipped with a 320x256 pixels InSb FPA sensor 
was used for the experiments. The IR camera has 
a maximum frame rate at full resolution of about 
100 Hz, while the integration time can be set 
between 10 pis and 20.000 pis, with steps of 1 pis. 
The integration time (r) used for acquiring the 
infrared images actually determines the maximum 
equivalent sampling frequency (f e , max ) that can be 
reached by the MTDA method, being: 



Eq.l 



The measurement principle of MTDA is 
schematically illustrated in Fig.l. 




Trigger out 
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Fig. 1: Schematic illustration of the MTDA method. 

Signals used for IR image acquisition are 
generated by means of three digital ports of an 
ATmega328 programmable digital micro- 
controller. The trigger output port of the IR 
camera (trigger OUT) generates a time base 
constituted by a train of impulsive signals having 
a period equal to the sampling time of the IR 
camera. When a rising edge occurs in the trigger 
output, a signal is generated by the 
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microcontroller and sent to the trigger input port 
(trigger IN) of the IR camera, starting the 
acquisition. A second signal is sent to the power 
device throughout a gate driver to start the heating 
event. These two signals are synchronized using 
the time base of the IR camera, thus providing an 
effective matching between the images 
acquisition and the transient event. 

To account for possible variation of the ambient 
temperature during the test (actually the number 
of cycles required for the reconstruction of the 
thermal transient can be quite high, depending on 
the dynamics of the event), a further single frame 
acquisition is triggered at the end of each cooling 
period and used for ambient temperature 
compensation. Hence, temperature increments 
over the initial temperature (delta_7) before the 
start of each heating cycle were considered. All 
the frames acquired at each trigger event are used 
at the end of the measurement procedure to 
reconstruct the time and space behavior of the 
temperature over the device. 

The timing sequence is schematically illustrated 
in Fig. 2, in which a) represents the switching 
time of the device under test, b) the timing of 
images acquisition and c) the trigger signals used 
for synchronization. 



delay (i.e., the equivalent sampling time), then the 
number of iterations will be: 




Fig. 2: MTDA for slow transient events (not in 
scale): a) switching time of the device under test, b) 

timing of IR images acquisition and c) trigger 
signals used for synchronization. D = heating time, 
R = repetition rate (cooling time), dt = time delay, T 
= pre-trigger used for ambient temperature 
compensation. 

If FR is the frame rate of the camera, at each 
iteration N infrared images will be acquired (note 
that, due to the time delay, sometimes one image 
lies outside the duration of the event and can't be 
utilized for the reconstruction): 



N = DxFR 



Eq.2 



where D represents the time duration of the 
transient event. The number of iterations depends 
on the time delay, and hence on the required 
equivalent sampling frequency. If f e is the 
equivalent sampling frequency and dt the time 



M 



fe _ 1 

FR dtxFR 



Eq.3 



Each iteration is shifted of dt allowing the 
complete reconstruction of the transient event that 
must be replicated M times. This holds as long as 
D is greater than the frame rate (i.e., for slow 
transient events). 

For fast transient event (i.e., when D is smaller 
than the frame rate of the IR camera) , the MTDA 
method can be still used by acquiring, at each 
iteration, only a single image instead of N images 
(Fig. 3). In this figure the shaded rectangle 
actually represents the IR camera integration time, 
which, according to (1), must be smaller than the 
time duration of the event. Of course, in this case 
the duration of the reconstruction process will be 
much longer. 



4^ 



Fig. 3: MTDA for fast transient events (not in 
scale): a) switching time of the device under test, b) 

timing of IR images acquisition and c) trigger 
signals used for synchronization. D = heating time, 
R = repetition rate (cooling time), dt = time delay, T 
= pre-trigger used for ambient temperature 
compensation, t = integration time. 



3. Validation of the proposed method 

To validate the MTDA method, a laser-cut 
planar resistor was chosen as reference target. The 
test circuit employed to carry out the experimental 
tests is schematically illustrated in Fig. 4a. The 
resistor (Fig. 4b) was heated cyclically from 
ambient temperature up to about 37 °C by means 
of pulsed currents of D = 1600 ms duration. No 
paint was applied onto the resistor to uniform the 
target area emissivity, which was estimated to be 
about 0.95 by separate tests. The repetition rate of 
the trigger pulse (R) has been chosen to be 1 min 
in order to allow the device to dissipate the excess 
temperature increase before a new trigger pulse is 
applied. A further single frame acquisition is 
triggered at the end of each cooling period and 
used for ambient temperature compensation, as 
already explained in the previous section. The 
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relatively slow dynamics of resistor's heating 
allows images acquisition to be performed also in 
normal operation mode (i.e., without time delay). 
All computations were carried out using Matlab™. 

Fig. 5 shows selected infrared images taken at 25 
Hz frame rate during one heating cycle (the 
complete sequence includes 40 images) .These 
images were taken without introducing any time 
delay (i.e., dt =0 ms) and served as reference for 
verifying the effectiveness of the proposed 
method and assessing the uncertainty of the 
reconstruction process. 

As expected, it can be observed that heat is 
generated suddenly as the current flows into the 
resistor, and eventually spread across the whole 
area of highly conductive material, with the 
hottest elliptically-shaped region located in the 
inner part of the resistive element. It can also be 
highlighted that the two conductive films used to 
drive current to the resistor element do not 
dissipate power significantly, because of their low 
resistance. In addition, this part of the resistor has 
a lower emissivity coefficient with respect to the 
black target area. 




Fig. 4: a) Test circuit for slow dynamics analysis 
and b) reference planar resistor used to carry out 
the experimental tests. 



1 


Q 




1 


\ 


■ 


y 


160 ms 


320 ms 




480 ms 




640 ms 


800 ms 


U 


B 


1 


y 


1 


B 




960 ms 


1120 ms 


1280 ms 


1440 ms 


1600 ms 



Fig. 5: Current-driven heating of the planar resistor 
(one cycle). 



The next figure (Fig. 6) reports normalized 
temperature profiles (delta_7 increments) of the 



hottest point of the resistor measured throughout 
five cycles, with no time delay. From this plot it 
can be realized that the repeatability of the 
observed event is quite good, being the average 
standard deviation (21 mK) comparable to the 
NETD value of the used IR camera. 
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Fig. 6: Normalized temperature profiles of hottest 
point of planar resistor (five cycles, dt = 0 ms). 



The next experimental tests were carried out by 
using the MTDA method to capture the current- 
driven heating of the resistor with an equivalent 
sampling rate of 200 Hz (dt = 5 ms) and 1000 Hz 
(dt = 1 ms), respectively. In the first case, the 
complete sequence includes 320 infrared images, 
while in the second one the number of infrared 
images captured is 1600. Each test was repeated 
five times and average heating profiles of selected 
points of the target area were computed. Fig. 7 
compares the reconstructed average heating 
profiles of the middle point (hot spot) of the target 
area with that obtained by using the infrared 
camera in normal operation mode (at 25 Hz frame 
rate). It can be observed a very good agreement 
among the three profiles. 




Time (ms) 

Fig. 7: Reconstructed heating profiles of hottest 
point of target area, resampled at 200 Hz (dt = 5 
ms) and 1000 Hz (dt = 1 ms). Values were measured 
by means of a cooled infrared camera operating at 
25 Hz frame rate. The circles show temperatures 
recorded in normal operation mode. 
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Fig. 8: Statistical deviation between average values 
of the heating profile of the hottest point of the 
target area sampled at 200 Hz and third order 

polynomial fitting of the heating profile measured 
in normal operation mode. 



Equivalent sampling rate: 1000 Hz 
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Fig. 9: Statistical deviation between average values 
of the heating profile of the hottest point of the 
target area sampled at 1 kHz and third order 
polynomial fitting of the heating profile measured 
in normal operation mode. 



In order to assess the measurement accuracy of 
the reconstruction process, the statistical residuals 
between the average values of the heating profiles 
obtained by MTDA and the third order 
polynomial fitting of the heating profile measured 
in normal operation mode were evaluated. Results 
are shown for an equivalent sampling rate of 200 
Hz and 1000 Hz in Figs. 8 and 9, respectively. An 
average error of (-0.08 ± 0.05) °C and (-0.06 ± 
0.08) °C was obtained, with no evidence of 
systematic effects related to the increasing of the 
equivalent sampling rate. 



4. Application to reliability 
assessment of IGBT power device 

The on-state voltage of power electronic devices 
tasked to perform several billions of repetitive 



switching operations during their life is subject to 
modifications caused by the progressive ageing of 
semiconductor and metal layers. Endurance tests 
are the traditional way to monitor these changes 
in order to assess device reliability. However, 
they are very time expensive, requiring even 
months of uninterrupted laboratory tests. An 
interesting alternative is to assess the reliability 
through a suitable mathematical model. 
Specifically, coupling the results of a very fast 
thermodynamic analysis with a reliability model 
based on the Coffin-Manson law, device 
degradation over time can be estimated and the 
level of reliability evaluated [10]. The proposed 
MTDA method was applied to investigate the 
temperature distribution over the source metal 
surface of an IGBT power device during a 
switching cycle. According to the test circuit of 
Fig. 10, in each switching cycle the device is 
turned on for 6 ms to generate a fast heating, and 
then turned off for 60 s, a time sufficient to allow 
the device to reach again the initial temperature. 
The time delay was set to 100 pis (corresponding 
to/e = 10 kHz), therefore it was necessary to carry 
out 60 iterations to completely reconstruct the 
heating event. A 1.6 mH inductor is exploited to 
limit the device current during the switching 
cycle. The device under test (DUT) was 
previously subjected to a suitable accelerated 
ageing procedure, then, the package molding was 
chemically removed to expose the source metal 
surface, as shown in Fig. 1 1 . 





Gate Driver 




1 





Fig. 10: Test circuit for IGBT fast dynamics 
analysis. 

Since the exposed surface encompasses regions 
made of different materials (aluminum, oxides, 
passivation etc.), a suitable emissivity map is 
needed to obtain the true temperature mapping. 
This was done by passively heating the device at 
constant temperatures. Moreover, in order to 
focus the investigation exclusively on the region 
of interest, captured images were post-processed 
by applying a mask exploiting the different 
emissivity of the target area with respect to the 
outer part of the device, which is made of silicon 
oxide. Before further computation, each infrared 
image was masked and corrected for the actual 
emissivity. No microscope lens was used for this 
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test, hence the spatial resolution is poor because 
of the small dimension of the device. 




Fig. 11: Live IR image of DUT 



Fig. 12 reports the reconstructed normalized 
temperature profiles measured in avalanche mode 
on one point of the DUT surface. Typically, a 
maximum deviation of about 1.1 °C was observed 
at the peak point. This value is comparable to the 
measurement accuracy of the IR camera {+1-1% 
or +/-1°C) as stated by the manufacturer. 
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Fig. 12: Normalized temperature profiles measured 
in avalanche mode on one point of the DUT surface 
(four cycles, dt = 100 us). 

Fig. 13 shows reconstructed temperature 
distribution maps of the DUT surface (in the 
figure, only images taken every 400 pis are 
shown) . 

Though the spatial resolution of the images is 
rather low because of the lens system (50 mm, 
F/2) used to perform these preliminary tests, 
nevertheless, the obtained results reveal the 
potentiality of MTDA method in mapping thermal 
transient events at full frame with outstanding 
time resolution. 

As far as the specific application is concerned, 
the ability to track the thermal evolution of IGBT 
power device allows important information to be 
obtained. For example, hot spots due to aging or 
failures as well as non-homogeneous current 



density distributions and thermal instabilities can 
be easily detected. 




Fig. 13: Reconstructed temperature distribution 
maps of the DUT surface, sampled at 10 kHz 
equivalent frame rate (dt= 100 jis). 



5. Conclusions 

A novel method able to extend the area of 
application of standard infrared systems to the 
investigation of very fast replicable 
heating/cooling transient events has been 
presented. The method basically relies on 
triggering multiple time-delayed acquisitions 
(MTDA) of the observed thermal phenomenon, 
by means of a very precise and stable 
programmable digital micro -controller and by 
reconstructing the time domain IR sequence using 
the frames acquired at each trigger event. As a 
result, the time resolution of full frame IR 
sequences is increased by two order of magnitude 
or more with respect to that of a standard infrared 
camera. Preliminary results obtained on a power 
electronics device with the molding package 
removed are very promising. The presented 
method, originally developed to assess the 
reliability of power devices, could find 
application in other areas as well. 
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Abstract 

Magnetic Resonance Imaging (MRI) thermometry 
is a non-invasive technique appropriate for 
temperature distribution monitoring in tissues. 
This technique, through a multidirectional scan 
without use of ionizing radiation, estimates the 
temperature variation thanks to the thermal 
dependence of several MRI parameters (i.e., Tj 
and T 2 relaxation times, proton resonance 
frequency shift and others). The use of T } shows 
some advantages, such as: low sensitivity to 
patient motion, high thermal sensitivity and 
linearity; on the other hand, the thermal 
sensitivity strongly depends on the organ tissue. 
The aim of this study is to assess the sensitivity of 
MRI thermometry using two Tj-weighted 
sequences (i.e., IRTF and SRTF) and an 1 .5-T MR 
scanner on healthy swine pancreases undergoing 
Laser Ablation (LA). The effect of temperature 
variation on the signal intensity was analyzed. 
The reference temperature was measured by MRI- 
compatible fiber optic sensors (fiber Bragg 
grating sensors). The sensitivity of the proposed 
techniques were estimated and compared. The 
IRTF sequence provides higher temperature 
sensitivity than the SRTF one (1.8 ±0.2 °C 7 vs 
1.4 ±0.1 °C~ 1 ) but the correlation between signal 
intensity and temperature variation is lower (R=- 
0.66 vs -0.97). Results show that the proposed 
technique may be adequate for temperature 
monitoring during LA. 

1. Introduction 

Laser Ablation (LA) is a minimally invasive 
procedure employed to treat tumor. Its safety and 
efficacy could be improved by performing a real- 
time temperature monitoring, hence several 
techniques have been investigated during the last 
thirty years. These techniques are usually divided 
in non-invasive approaches and invasive ones [1]. 
The interest on non-invasive approaches is 
motivated by the advantages related to the non- 
invasiveness and to the ability of these techniques 



to provide a tridimensional map of tissue 
temperature. The most investigated non-invasive 
approaches are: 1) MRI thermometry, which is 
based on the sensitivity of several MRI 
parameters on temperature [2]; 2) CT -based 
thermometry, which is based on the influence of 
temperature on images obtained by Computed 
Tomography scans [3,4]; and 3) ultrasound-based 
thermometry, which is based on the dependence 
of several ultrasound parameters on temperature 
[5]. 

MRI thermometry shows some advantages, such 
as high thermal sensitivity, low sensitivity to 
motion, good linearity and it does not employ 
ionizing radiations [2] . 

The first investigation about the influence of 
temperature on MR parameters was conducted by 
Bloembergen et al. [6] in 1948, and Jolesz and 
coauthors developed the technique of MRI-guided 
LA [7]. Since the study of Jolesz (1988), a big 
research effort has been dedicated to assess the 
feasibility of MRI thermometry for temperature 
monitoring during LA, and among the several MR 
parameters, T t relaxation time and proton 
resonance frequency (PRF) have been considered 
the most attractive [8] . 

The aim of this work is twofold: the sensitivity 
analysis of MRI thermometry for temperature 
monitoring of pancreatic tissues during LA; the 
comparison between the characteristics of MRI 
thermometry using two T r weighted sequences. 
The reference temperature for MRI thermometry 
characterization has been measured by MRI- 
compatible temperature sensors (i.e., fiber Bragg 
grating, FBG, sensors). To the best of our 
knowledge, this study represents the first 
investigation of MRI thermometry on pancreatic 
tissue, and this is important because of the 
influence of the tissue histological characteristics 
on Ti thermal sensitivity. A further novelty of this 
work is related to the use of new settings for the 
two employed T r based sequences. 
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2. Theoretical background 

A number of MR parameters depend on 
temperature [2]: proton density, T t and T 2 
relaxation time, magnetization transfer and proton 
resonance frequency (PRF) among others. 

As is well known, after the transmission and the 
successive interruption of a radio frequency signal 
to the medium, proton spin reacquires the initial 
position and direction. The mechanisms of 
relaxation of the spin proton to recover the 
equilibrium status are described by two temporal 
parameters, i.e., spin-lattice relaxation time Ti and 
the transversal relaxation time T 2 . Therefore, they 
provide information about the structure of the 
tissue. For application in temperature mapping the 
time T 2 is not attractive because its thermal 
sensitivity is lower than the one of Ti [9] and 
because it can be masked by other factors [10]. 
Similarly, variations in the proton density can be 
influenced not only by temperature but also by 
other parameters, such as Therefore, mostly 
used parameters are Ti and PRF, because of their 
good sensitivity to temperature variation. The 
main limitation of the first parameter (i.e., Ti) is 
the sensitivity dependence on tissue; on the other 
hand, the negligible sensitivity to movement is its 
main advantage. PRF is not tissue-dependent and 
has a linear dependence from temperature in 
ranges wider than those of Ti; however, PRF 
shows high sensitivity to movement, which can be 
an issue during in vivo applications [2] . 

We report in detail only the principles of T r 
based MRI thermometry because we used T r 
weighted sequences during experiments, as 
described in Section 3 . 

Ti of water protons is a measurement of the time 
required to protons to return to the initial 
equilibrium conditions, after the interruption of 
the radio frequency signal, thanks to the transfer 
of energy to the surrounding micro -environment. 
Spin relaxation results from the interaction of the 
dipoles with the lattice, which is due to their 
translational and rotational motion. The 
temperature (T) dependence of this motion causes 
the Ti increase with T, which can be expressed as 
follows [2]: 

T^e^-Wj Eq.l 

where E a (Ti) is the activation energy of the 
relaxation process, k is the Boltzmann constant, 
and T is the absolute temperature. 

This relationship can be linearized in a wide 
temperature range (e.g., from 30 °C up to 70 °C): 

Tj-T lref =m\T-T ref ) Eq.2 

where T ref is the reference temperature, m is the 
thermal coefficient, empirically determined for 



each tissue [11], and T lref is the T x value at a 
temperature T ref . Non-linear effects can occur if 
the tissue properties change, for instance they can 
change because of the tissue coagulation during 
LA. Furthermore, the thermal sensitivity of T t 
depends on the tissue type (e.g., 1-2% -°C A in 
liver, 1.4%-°C _1 in bovine muscle, and 0.97%-°C _1 
in fat), therefore the previous knowledge of the 
thermal coefficient of each tissue is essential to 
obtain a temperature map. The research group 
headed by TJ Vogl is performing a deep analysis 
of T\ weighted images during LA of liver 
metastasis, and of liver-mimicking acrylamide gel 
phantom [12,13]. 

In order to monitor the tissue temperature during 
LA and during MRI scans, we employed fiber 
Bragg grating (FBG) sensors, which do not 
present artifacts because of their immunity from 
electromagnetic interferences and their MRI- 
compatibility [14]. They consist of a periodic 
variation in the reflective index of the fiber core. 
The grating is characterized by its spatial period, 
A, and the effective reflective index, n eff . When 
interrogated with a broadband radiation, a narrow 
range of wavelengths is reflected, the other ones 
are transmitted. The reflected wavelengths range 
is centered around a specific value, i.e., the Bragg 
wavelength (X B ), expressed as follows: 

A B =2-A-rj eff Eq.3 

The X B variation (AA, B ) depends on the variation 
of temperature, AT, and on mechanical strain, 8. 
During LA, considering negligible the strain, the 
AX B can be considered only dependent on the 
temperature increments and therefore it provides 
an indirect measurement of temperature increase 
during the treatment. 

3. Experimental setup, data 
acquisition and processing 

LA was performed on swine pancreases by a Nd: 
YAG laser (X=1064 nm) that conveyed the 
radiation into a fiber applicator. Each organ was 
placed into an MRI-compatible box (10x10x2 cm) 
composed of (poly)dimethylsiloxane (PDMS), 
necessary to insert the applicator and FBGs and to 
control their relative distances during LA. Figure 
la shows a 3-D view of the box and the positions 
of the laser applicator and the four fibers Fl, F2, 
F3, F4, each of them equipped with 3 FBG 
sensors. The FBG sensors have a sensing element 
with 1 mm of length. The laser applicator was 
inserted into the tissue at a depth of 2 mm and at 
the center of the box. 

Three FBGs (i.e., Gl, G2 and G3) housed in 
each optical fiber, are placed at a distance of 2 
mm each other; considering that the x-y plane of 
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the applicator is perpendicular to the plane of the 
optical fibers, the central FBG -i.e., G2- is placed 
at the same quote of the applicator (fig. lb). The 
four fibers were placed at the same four positions 
in all the tests; these positions correspond to the 
center of the four Region of Interests (ROIs) 
which were selected to analyze the images. 




Fig. 1: a) polymeric box used to control distances 
between laser applicator and optical fibers Fl, F2, 
F3, F4; b) arrangement of FBGs (Gl, G2, G3) inside 
each fiber; c) pancreatic tissue inside the box placed 
within MR coils. 



The heating process was monitored by an 1 .5-T 
MR scanner (Siemens Magnetom Avanto, 1.5 T), 
using two T r weighted sequences: inversion 
recovery turbo flash (IRTF) and saturation 
recovery turbo flash (SRTF), fig. lc. The 
variation of T 1? due to the increase of temperature 
(Eq. 2), causes the changes of image signal 
intensity, AS, in T r weighted images. In order to 
obtain the thermal dependence of AS, images 
were analysed in ROIs centered in 
correspondence of the FBG sensors. The FBGs 
output (Bragg wavelength, X B ) has been recorded 
by an optical spectrum analyzer, OSA (Optical 
Sensing Interrogator, sml25, Micron Optics). The 
OSA records the wavelength shift due to the 
temperature variation. The temperature increment 
during laser ablation is estimated by the 
calibration curve of the FBG sensors, which were 
calibrated in a wide range of measurement (i.e., 
from 18 °C up to 100 °C). 



4. Results and Discussion 

In order to obtain the relationship between AS, 
and the temperature increase, AT, the averaged 
intensity of each ROI was synchronized and 



correlated with the temperature values monitored 
by the FBGs. 

First of all the noise of the images obtained 
using the two sequences has been analysed. We 
considered the standard deviation, std, of the 
pixels intensity in a ROI with no nuclear MR 
signal. The std was calculated considering a ROI 
with about 5000 pixels in two images, the first 
one obtained using SRTF and the second one 
using IRTF. The STRF image showed lower noise 
(std=0.67) than the IRTF image (std=l .2). 

Regarding the analysis of pixel intensity (S) in 
the image representing the tissue, it must be 
considered that the std increases because of the 
tissue dishomogeneity. A typical histogram 
obtained considering a ROI of about 5000 pixels 
representing a pancreas is shown in fig. 2 for the 
two sequences. 




90 95 100 105 0 110 115 120 



Fig. 2: histograms of the signal intensity of a ROI 
representing a pancreas obtained using the two 
sequences. 

The std values of the image obtained by SRTF is 
slightly lower than the one obtained by IRTF (i.e., 
4.8 and 5.1 for SRTF and IRTF, respectively). 
The tissue dishomogeneity should influence the 
performances of the MRI-thermometry. 

After this preliminary analysis has been 
investigate the changes of signal intensity during 
the treatment, caused by the temperature increase. 
In particular, the value of AS, which is the ROI 
value difference between the current and the 
reference image acquired at the beginning of the 
trial, has been calculated for all the MRI scans 
performed during the procedure. Figure 3 shows 
AS as a function of AT. 
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Fig. 3: Signal intensity vs temperature for the two 
sequences. The best fitting lines are also shown. 

Figure 3 shows that the use of IRTF allows 
increasing the sensitivity: the slope of the best 
fitting line, which represents the thermal 
sensitivity, is 1.8±0.2 °C 1 using IRTF and 
1.4±0.1 °C 1 using SRTF. On the other hand the 
correlation with AT improves using SRTF, as 
shown in fig. 3 and confirmed by the values of 
correlation coefficient, R (i.e., -0.97 and -0.66 for 
SRTF and IRTF, respectively). 

The AT has been calculated by the calibration 
curve of the FBG sensors. All the FBG sensors 
were calibrated using an experimental set up 
composed of: i) an optical spectrum analyzer 
(Bragg Fiber Sensing, FS2200 8 CH, Sequoia 
Technology Group Ltd), which records the output 
(Xb) of the FBGs, it) a PC to record the output 
data from the analyzer, Hi) an oven to control the 
temperature, iv) a module to acquire and record 
the signal provided by four K-type thermocouples 
(FX106-4-2, YOKOGAWA®), v) a PC to process 
data. 



The calibration was performed between 18 °C 
and 100 °C. The mean value of the measures 
obtained by the four thermocouples was used as 
reference for the FBGs calibration. The change of 
the Bragg wavelength (AX B ) for each sensor 
during the whole process was recorded and 
synchronized with the temperature provided by 
the thermocouples in order to obtain the 
calibration curves of the FBGs by fitting the AX B 
with T. 

The curves have been obtained by a linear 
regression of the data. A typical curve was: 



AX n = 0.01 0-T- 025 



Eq. 



Therefore the sensitivity of the FBGs was about 
0.010 nm- 0 C _1 . 

Figure 4 shows a comparison between the 
temperature increases measured by FBGs (T ref ) 
and the ones estimated by the two sequences 
(Tirtf and T S rtf) obtained by considering the two 
best fitting lines reported in fig. 3. 



9 10 11 12 13 14 15 



| AT r- 



0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 

Fig. 4: comparison between the temperature 
increases measured by FBG and by the two 
sequences. For both the sequences data obtained 
during the LA performed on two pancreases are 
reported. 



Data reported in Fig. 4 show that for the 14 
scans performed during LA of two pancreases the 
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differences between the AT estimated by SRTF 
sequence and the AT estimated by FBGs are 
lower than the ones between the AT estimated by 
IRTF sequence and the AT estimated by FBGs. 
Considering the high accuracy of FBG 
measurement (i.e., -0.1 °C), the T estimations 
performed by SRTF sequence have a smaller 
measurement error than the ones measured by 
IRTF. This consideration is also confirmed by the 
Root Mean Squared Error for the two sequences 
calculated as: 



RMSE IRTF =^ 



N , 



AT, 



ref 



N 



Eq. 5 



and 



RMSE SRTF =\ 



Y.{AT ref -AT SRTF } 



N 



Eq.6 



In fact the value of RMSE IRTF resulted bigger 
than the one of RMSE S r T f (i.e., 8.3 °C vs 1.8 °C) 

The improvement obtained by performing the 
estimation by SRTF has also confirmed by the 
Bland Altman analysis: the narrow range defined 
by the limit of agreements (LOAs) for SRTF (i.e., 
from -2.3 °C to 2.1 °C) confirms the better 
precision of these measures than the ones 
obtained by IRTF (i.e., from -9.4 °C to 9.1 °C); 
also the mean of differences (MODs) show a low 
value (i.e., 0.1 °C for both the sequences). 

5. Conclusion 

Summing up, the sensitivity during LA on 
healthy porcine pancreas of two T r weighted 
sequences has been evaluated by comparing the 
AS with the AT measured by MRi compatible 
sensors based on fiber Bragg technology. The 
SRTF sequence has shown lower sensitivity than 
the IRTF one; on the other hand it shows a better 
precision than IRTF. 

In conclusion, we have investigated the 
feasibility of MR-based thermometry for 
monitoring pancreatic tissue temperature during 
LA, and the two Ti -weighted sequences can be 
assessed positively considering the MOD and 
LOA values especially for the pancreas SRTF 
sequence: MOD SF rt = -0.1 °C and LOA SF rt =[- 
2.3; 2.1] °C. Moreover, SRTF and IRTF are 
standard sequences, therefore commonly used in 
clinical practice for the traditional analysis. 

Real time MR-based thermometry could 
represent a valuable tool to lead the physician in 



the optimal laser settings during thermal ablation 
procedure. This technique is particularly attractive 
for LA due to the laser MR-compatibility. 
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Abstract 

Rolling contact fatigue (RCF) plays a critical role 
in railway components, and the characterization of 
materials used, in terms of RCF life, is still an open 
task, made complex by the interactions of different 
phenomena, such as wear, cyclic contact, and 
ratcheting. 

The presented case study regards a bi-disk test 
bench, used to evaluate the RCF behaviour of 
wheel and rail steels. In the test, a wheel steel 
specimen rotates against a rail steel specimen, 
while pressed, the one against the other, by a 
constant force . 

A numerical simulation, based on a multiple 
source damage model and on parameters obtained 
from direct measurements during the test, is then 
used to determine properties of the material. 

The contact surface, in particular, plays a key 
role in determining the evolution of an RCF life 
test, since it has a direct impact on the pressure 
exerted and can change during the test, due to 
wear. 

The procedure proposed aims at using vibrations 
of a test bench during RCF -life tests to identify the 
contact topology, specifically understanding when 
damage phenomena cooperate to cause a quick 
flattening of the surface, and when this process is 
complete. 

1. Introduction 

Characterization of RCF life of materials is critical 
in railway applications due to the fact that cyclic, 
localized stresses lead to rolling contact fatigue 
(RCF) phenomena, which can cause severe 
damage and sudden failure of a component. 
Wear, also, is a common phenomenon occurring on 
both rail and wheels, but its interaction with RCF 
is not always detrimental: in specific conditions 
(uniformity, low rate, stability), it could optimize 
contact geometry, increasing RCF life, while, in 
non-optimal conditions, it could reduce RCF life. 
Moreover, the high shear stress has to be taken into 
account, leading to a complex combination of 
concurring damage phenomena [1,2]. Models 
describing cyclic plasticity are available in 



literature, although their results are not reliable if 
material parameters are obtained by laboratory 
tests in standard stress condition [3,4], especially 
due to a different propagation mechanism [5] . 
A reliable procedure to analyse RCF-wear 
interaction would be to carry on accelerated RCF 
tests on controlled slip ratio conditions with 
multiple specimens, interrupting the test at 
different times and performing destructive analysis 
on the specimen, which is costly and time 
consuming. 

The alternative proposed it to perform a single 
RCF test, monitoring different mechanical 
quantities to detect changes in the process. 

2. Materials 

The proposed procedure was applied to a pair of 
wheel and rail materials (ER8 EN13262 an UIC 
900 A respectively), as a case study. 

The wheel specimens were 60-mm wide and 10- 
mm thick cylindrical disks, while the rail 
specimens were 59.5 mm in diameter, a 10 mm 
thickness and had a crowning radius of 200 mm, to 
prevent any border effect on the contact surface. 

The mechanical properties of wheel rims and rail 
heads of the manufacturer are reported in Table 1 . 



Steel 


ER8 


900A 


Ultimate tensile stress 
[MPa] 


940 


930 


Monotonic yield stress 
[MPa] 


590 


470 


Cyclic yield stress 
[MPa] 


470 


390 


Necking [%] 


54 


26 


Elongation [%] 


17 


14 


Brinell hardness [HB] 


230-255 


296 



Table 1 . Mechanical properties of the steel used 



Three different contact pressure levels were 
tested, in two different sliding ratio conditions, for 
a total of 6 different rolling/sliding parameter 
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combination, here reported in Table 2. Vibration 
recordings from those with highest sliding ratio 
(S=3%) and with the highest nominal pressure 
(P=1500 MPa) were also examined using the 
proposed procedure. 





L 


R 


S 


P 


ID 


Contact 


Rolling 


Sliding 


Nominal 


load 


speed 


ratio 


contact 
pressure 




[N] 


[r.p.m.] 


[%] 


[MPa] 


A 


1,500 


516.0 (R) 
492.5 (W) 


a 
J 


1 ,1UU 


B 


2,490 


516.0 (R) 
492.5 (W) 


3 


1,300 


C 


3,830 


516.0 (R) 
492.5 (W) 


3 


1,500 


D 


3,830 


511.0 (R) 
497.5 (W) 


1 


1,500 


E 


2,490 


511.0 (R) 
497.5 (W) 


1 


1,300 


F 


1,500 


511.0 (R) 
497.5 (W) 


1 


1,100 



Table 2. Test parameters (W=wheel, R=rail) 



3. Methods 

3.1. Measurement system 

The test bench used is a bi-disk machine 
dedicated to study interactions between two 
components subjected to cyclic contact in different 
load conditions, already presented in previous 
works by the authors [6] , and whose general layout 
could be found in Figure 1 . 

Specimen 2 Mobile mandrel 







ffl. 
1 1 


1 w 

♦ I * i\ 




Fixed mandrel / 
Specimen 1 






Sci 


vo - hydraulic actuator 



Figure 1: Test bench layout 

The contact load, up to 70 kN, is maintained 
thanks to a servo -hydraulic actuator enabling the 
sliding of one of the mandrels, while two 
independent 33 kW engines provide the specimens 
rotation. Both engines and the actuator are 
controlled to perform tests at given slip ratio, 
rotating speed, contact force and engine torque. 



On both mandrel supports, piezo-accelerometers 
were positioned, one in the vertical and one in the 
horizontal plane, both normal to the rotation axis 
and having a sensitivity of 0.98 V/(m/s 2 ) and a 
linear bandwidth in the 5-20 kHz range. 

The signals from the two accelerometers, as well 
as those from a torque sensor, positioned on the 
sliding mandrel, were acquired by means of a 
configurable data acquisition system at a 5 kHz 
synchronous sampling frequency. 




Figure 2: specimens mounted on mandrels. The 
accelerometers measuring axes are highlighted 

An L101b-2k BASLER linear array camera, in a 
setup leading to a resolution of 1 //m on the 
specimen mid surface, in sync with the bench 
rotation, was also used to monitor contact surface 
conditions of both specimens, but, due to 
lubrication, it was only used at regular intervals 
while the test was on halt. 

Every 2xl0 5 cycles the test was halted to perform 
weighing with a precision scale (0.01 g resolution), 
to evaluate weight loss due to wear. 

At the end of each test samples were cut in mid- 
thickness, polished and etched (using 2% Nital) 
and inspected using an optical microscope to 
measure the plastic flow. 

Following this, RCF crack paths were also sought 
using a LEO EVO-40XVP scanning electron 
microscope with a microprobe, and HVio micro- 
hardness was measured at different depths on the 
cross-section of each sample. 

3.2. Mechanical Model 

While vibration analysis is commonly used to 
detect faults on gears and bearing components [7- 
10], its use in material characterization is limited, 
as this is usually carried on with standardized 
methods, generally involving destructive testing or 
evaluation at the end of test. 

In previous works [11] a damage indicator, 
estimated from vibrations and correlated to surface 
cracks was proposed. Following the same approach 
the two specimen interactions are synthesized, as it 
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is common in a modal analysis, using frequency 
transfer functions (FRF) between torque T and 
accelerations in the x and y directions, as depicted 
in Figure 3 . 

The idea behind this approach is that any 
variation in the material or geometry of the 
samples would be reflected by a simultaneous 
change in the FRFs values, helping identifying 
when crowing gets flattened and the process 
duration. 




Figure 3: Simplified FRF based model 



FRFs were numerically computed as the Ha/b 
cross-spectrum between a signal a and a signal b, 
divided by the autospectrum of signal b. In 
particular, the quantities monitored were the 
transmissibility H y / X , between the y acceleration on 
the fixed mandrel support along the vertical axis 
and the x acceleration of the moving mandrel 
support along the pneumatic actuator axis, as well 
as the rotating inertances Hx/t and H y / T between the 
x and y accelerations and the torque T measured on 
the rotating axis of the sliding mandrel. To provide 
a complete set of transfer function indicators, the 
reciprocals to those three quantities were also 
computed, resulting in the computation of H x / y , 
H T /x and H T / y , 

Though all these quantities are related to 
properties of the system, such as masses and 
stiffnesses, their connection to any proper 
mechanical parameter would require an overly 
complex multiple degree of freedom model, which 
would exceed the scope of this research. Therefore 
their absolute values are not considered 
meaningful per se, and only their variation from 
their initial value of each test has been taken into 
account. 

3.3. Preliminary tests 

Preliminary tests were carried on to select which 
quantity to monitor and the frequency range in 
which the changes sought have a higher impact. 

To avoid misreadings due to external or 
uncorrelated sources of vibrations, the test bench 
was monitored while running in jog mode, as well 
as in full operational mode, but without any load 



applied between the samples. A power spectral 
density analysis of x, y and T signals, recorded for 
about 100 s, and averaged using a 1 Hz spectral 
resolution, revealed that their frequency content 
was located below 400 Hz, while the same analysis 
performed in the final phases of standardized RCF 
tests revealed vibrational phenomena up to 1 ,000 
Hz. 

An experimental modal analysis of the pneumatic 
actuator, also, pointed out a resonance frequency 
lower than 100 Hz, therefore, the spectral 
quantities monitored were treated using a band- 
pass filter between 400 and 1,000 Hz. 

A comparison between the quantities monitored 
during a full test of an unrelated rail/ wheel couple 
[11], here reported in Figure 4, pointed out how the 
FRF between torque and acceleration along the 
linear actuator axis, and its reciprocal, displayed a 
sensible change, coherent with other damage 
indicators, and a behaviour less affected by spikes 
associated with accidental impacts or other 
impulsive events due to the uncontrolled 
environment in which the bench operates. 
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Figure 4: comparison between H y / X , Ht/ x and HT/y, 
other indicators pointed out a wear rate change 
after 60,000 and 100,000 cycles [11] 

3.4. Vibration analysis 

Using the acquired data, the transfer function 
Ht/ x (co) between torque and acceleration was 
computed averaging 5 windows of 1 s, thus 
producing a spectrum of 1 000 lines every 5 s . 

The initial modal properties of the system made 
up by the test bench and the intact specimens were 
assessed by averaging the Hx/ y (o3) values of the first 
2500 cycles, thus assumed as a set of reference 
values Hx/ x ref (co). 

To point out changes in the system, the difference 
dT/x(a>) between the current status and the reference 
value, for each frequency resolution ©i, has been 
computed, then the synthetic index D t has been 
evaluated as the root mean square (RMS) 
summation for all frequencies, as shown in 
Equation 1. 
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To further smooth the resulting signal, given the 
relative slow process, the RMS value was averaged 
on a 30 s window. The same procedure, here 
illustrated in Figure 5, was performed for the 
reciprocal Hx/ T , leading to an indicator D x / T , which 
was also monitored during the RCF tests. 



Torque and Acceleration Measurements 



5 s averaged FRF 



30 s averaged FRF 



400-1000 Hz Band pass filter 



Current FRF 



Reference FRF 



sliding ratio level (3%) report a steep increase in 
value (from 10 to more than 100 kg^m" 1 ), after few 
thousand cycles, followed by an abrupt decrease 
after a variable amount of cycles, then settling 
around a stationary level (30 kg^m" 1 ). 

The last test, with a 1% sliding ratio, displayed 
only a slight increase after 8,000 cycles to a 
stationary level (20 kg^m" 1 ). 

Such behaviour could be used to identify a time 
frame, characterized by higher D y /x levels, during 
which the damage process is clearly different than 
the initial or stationary stages: the approximate 
onset and duration of these periods are reported in 
Table 1. 



Difference between spectra 



30 s averaged RMS 



Figure 5: Signal processing procedure's step 

4. Results 



ID 


Onset 


Offset 


Duration 


[cycles] 


[cycles] 


[cycles] 


[s] 


A 


700 


4,700 


4,000 


480 


B 


3,800 


8,500 


4,700 


560 


C 


6,000 


12,700 


6,700 


800 


D 


-n/a- 


8,200 


-n/a- 


-n/a- 



Table 1: onset and duration of high Dy/T levels 



4.1. Vibrations 

The values of Dx/t for the four tests under scrutiny 
are reported in Figure 6. All three tests with a high 
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Figure 6: Dx/t indicators for tests A (black), B (red), C (green) and D (cyan) 



4.2. Weight loss 

The wear rate, in terms of weight loss as a 
function of cycle number, with varying nominal 
contact pressure P, is shown in Figure 7 and Figure 
8 for sliding ratios s = 1% and s = 3% respectively. 
Given the substantially linear relationship between 
weight loss and cycle number in the steady state 
period of the wear curves, the wear rates were 
calculated from the linear best fit of the 



experimental points. Their values increased as 
nominal contact pressure and sliding ratio 
increased. The wear rate in general was higher for 
the rail steel samples, despite the surface damage 
was more severe for the wheel steel ones: this 
means that wear mitigates the effect of surface 
pitting by removing damaged layers. This linear 
wear rate behavior suggests a stationary contact 
geometry already achieved before the first 
measurement (after 200,000 cycles), while the 
difference between 3% and 1% sliding ratio tests 
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points out a decreased wear rate associated with 
lower sliding condition. 
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the repeated on 20 images per sample and then the 

ratio between damaged and visible area was 
averaged and multiplied by the nominal thickness 
of the specimen. When no clear traces of the 
original surface of the wheel specimen was 
identified in the image, the whole surface (10.0 
mm) was considered as damaged. 
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Figure 7: 1% sliding test weight loss 
for rail (upper) and wheel (lower) 
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Figure 9: wheel contact surfaces before and after 
the first 200,000 cycles for test A 
Area worn by contact is highlighted. 
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Figure 8: 3% sliding test weight loss 
for rail (upper) and wheel (lower) 

4.3. Surface images 

Shape changes due to wear and RCF were 
assessed by measuring the worn portion of the 
wheel specimen as shown on images recorded by 
the linear camera before test and after 200,000 
cycles, here displayed in Figure 9 toFigure 12. 

The worn portion was identified by manually 
selecting the image area of the specimen surface 
where linear patterns, due to the specimen 
machining during production, were not present, or 
were only partially visible. The process was 
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Figure 10: wheel contact surfaces before and after 
the first 150,000 cycles for test B 
Area worn by contact is highlighted. 
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Figure 11: Wheel contact surfaces before and after 
the first 200,000 cycles for test C 
The area worn by contact is highlighted 




Figure 12: Wheel contact surfaces before and after 
the first 200,000 cycles for test D 
The area worn by contact is highlighted. 

4.4. Scanning Electron Microscope 

The inspection of specimen sections pointed out 
that the main damage mechanism, in both rail and 
wheel steels sample, was ratcheting, initiating 
cracks which propagate to cause RCF failure. 

Figure 13 shows the RCF crack path for various 
test conditions: the rail steel samples (a and c) were 



always less damaged than the wheel steel ones (b 
and d) in all tested conditions, and presented only 
surface cracks. In addition, the wheel steel sample 
also presented subsurface cracks. 

All surface and subsurface cracks analyzed 
display a shallow angle to the surface and follow 
the plastically deformed material during their 
growth. 




Figure 13: RCF crack path in: a) rail steel sample 
tested at s = 1% and P = 1,500 MPa; b) wheel steel 
sample tested at s = 1% and P = 1,500 MPa; c) rail 
steel sample tested at s = 3% and P = 1,300 MPa; d) 



wheel steel sample tested at s = 3% and 
P = 1,300 MPa 

4.5. Micro-hardness 

After tests were completed, micro-hardness 
measurements were taken at different depths for 
each specimen. The resulting hardness profiles are 
shown in Figure 14. 
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Figure 14: Micro -hardness profiles on the cross- 
section of rail and wheel steel samples 
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A higher hardening was observed in tests with 
higher sliding ratios and contact pressures, 
involving layers between 0.4 to 0.8 mm. 

5. Discussion 

The plastic behaviour was numerically simulated 
using a formulation by Mazzu for rolling contact 
[13], to predict plasticization in a two-dimensional 
plane-strain half-space subjected to contact 
pressure and friction, also taking into account wear 
as a concurring phenomenon removing layers from 
the surface. 

The micro-hardness profiles pointed out a 
behaviour associated with isotropic hardening, 
which allowed for a model simplification, 
supposing a linear relationship between yield stress 
and hardness increase [14]. 

The images taken after the first 200,000 cycles, 
displayed in Figure 9 to Figure 12, allowed to 
prove that wear flattened the rail specimen, thus 
increasing the contact surface and reducing the 
actual contact pressure exerted. 

Vibrations analysis, in particular of the first 
20,000 cycles reported in Figure 6, identified a 
timeframe in which the flattening process was 
occurring. 

Time and extent of the flattening allowed to 
replace the standard point Hertz contact in the 
model with a line contact surface, with a track 
width of 10 mm (8 mm for the 1% sliding ratio), 
occurring at the very beginning of the RCF 
simulations, as shown in Table 2. 



layers from the surface, preventing deep crack 
propagation. 
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[mm] 


[mm] 


[cycles] 


A 


1.785 


10 


4,700 


B 


2.113 


10 


8,500 


C 


2.440 


10 


12,700 


D 


2.440 


8 


8,200 



Table 2: onset and duration of high Dy/T levels 

The material properties needed for the simulation 
were calibrated by iteratively repeating simulation 
until the predicted displacements of material points 
due to plastic flow were in agreement with the 
experimental ones. 

The superposition of the calculated displacements 
on the observed grain boundaries, for two different 
test conditions, is presented in Figure 15. A good 
agreement between the two is evident. 
The simulation correctly predicted that the cracks 
and plasticized layers cannot be deeper than 40- 
50 jum, even though the ratchetting never stops, 
because wear removes plasticized and cracked 




Figure 15: Superposition of calculated and observed 

plastic strain bands in tests with s = 1 % and 
P = 1,100 MPa (left), and s = 1% and P = 1,300 MPa 
(right), after 2,000,000 cycles 

6. Conclusions 

A procedure to assess the presence, onset and 
duration of collaborative RCF and wear 
phenomena leading to contact shape alteration was 
proposed. The approach was applied to a wheel and 
rail steels characterisation procedure based both on 
rolling-sliding contact tests, on disk specimens, 
and numerical simulation. 

The advantage of the proposed method rests in its 
non-destructive and in-line nature, therefore 
allowing for a continuous assessment without 
halting tests or interfering with the specimens. 

This procedure completes the described full RCF 
failure life test, and allows to correct numerical 
models required to characterize the materials. It 
points out phenomena whose duration is too short 
to be taken into account by the same numerical 
model describing, for the whole RCF life of the 
specimen, the different damage phenomena and 
their interactions in a complex case as cyclic 
contact. 
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Abstract 

This paper describes the results of preliminary 
Monte Carlo simulations to investigate feasibility, 
performances and limitations of an original 
application of cosmic ray tracking to monitor 
static stability of civil structures, in particular 
historical buildings, where conservation 
constraints are severe and the time evolution of 
eventual deformation phenomena may be of the 
order of months or years. The simulated system 
consists of three or more particle detectors 
vertically aligned and linked together with the 
considered building. When a cosmic ray crosses 
all the measuring detectors, its track can be 
reconstructed and the reciprocal position of 
detectors can be determined. Continuous 
measurements allow to ascertain possible 
variations in detector positions related to 
deformations of the studied building. 

1. Introduction 

At the sea level, cosmic ray radiation is mostly 
composed of high energy muons that continuously 
hit the Earth's surface at a rate of 
10000 muon/(m 2 -min). Their energy distribution 
has a mean value of 3 GeV with a long tail at high 
energies and their angular distribution has its 
maximum value around the zenith direction [1] . 

Cosmic rays are largely exploited in nuclear and 
elementary particle physics for detector testing 
and calibration, and for the alignment of detectors 
in the very complex apparatuses used in this field 
[2]. Recently, cosmic ray muons are being 
increasingly considered as a source of natural, 
free, ubiquitous, high-penetration radiation of 
potential use for a number of applications: from 
monitoring in geology [3, 4, 5] to iron and steel 
industry [6], from security purposes [7, 8, 9, 10, 
11, 12] to monitoring alignment and stability of 
large mechanical and civil structures [13, 14, 15]. 

The present study regards preliminary Monte 
Carlo simulations to investigate the feasibility, 



performances and possible limitations of a system 
based on cosmic ray muon tracking, that monitors 
stability of the "Palazzo della Loggia", shown in 
Fig. 1, an historical civil building, seat of the 
municipal hall in the town of Brescia, in Italy. The 
top beam of its vaulted wooden roof has 
undergone a progressive deflection of 1 mm per 
year [16]. In particular, the standard measurement 
uncertainty, as a function of time, has been 
assessed, and the results have been discussed 
in comparison with a conventional monitoring 
system. Moreover, possible developments of 
the proposed technique have been considered 
and discussed, in order to improve the 
performances of the measurement system. 




Fig. 1: A picture of the Palazzo della Loggia, located 
in Brescia, Italy. 



2. The case study: the Palazzo 
della Loggia 

2.1. Overview of the case study 

The Palazzo della Loggia in the city of Brescia, 
Italy, presently hosting the Municipal Hall, was 
completed in 1574. Since then, it has cumulated a 
long sequence of injuries, transformations, 
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repairing interventions, some of which have 
generated considerable problems of structural 
stability of the building. The wooden vaulted roof, 
in particular, was completely reconstructed in 
1914: its maximum elevation is 16 m, and its 
shape is that of an upside-down ship, with planar 
rectangular sides of about 25 m and 50 m 
respectively. 

The structural architecture of the vault consists 
of principal truss wooden arches and simple 
secondary arches; both are connected at the top by 
a truss-made wooden beam. Immediately after its 
construction, the structure underwent progressive 
deformation of the longitudinal top beam and of 
the key points of the connected arches. In 
particular, the progressive deflection of the top 
beam was measured to be 190 mm in 1923, 
520 mm in 1945, 800 mm in 1980 [16]. 

2.2. The traditional monitoring 
approach 

Keeping the wooden vaulted roof of the Palazzo 
delta Loggia under control consists in measuring 
the possible displacements of several constituent 
elements with reference to a coordinate system 
fixed with the masonry structure of the building. 
The measurements should be taken at a frequency 
of days/weeks, for a total duration of several 
years. 

The systematic campaign of investigation by 
traditional measurement systems has been 
continuously performed for more than ten years, 
from 1990 to 2001 [16]. On four, out of the seven 
principal arches, three couples of wires, 2.0 mm 
in diameter, one made of ordinary steel and the 
other made of invar, were stretched between 
symmetric points at three different levels: Al - 
Bl, at the point of connection of the arches with 
the building structure; A2 - B2 and A3 - B 3, on 
the arch reins. These points are shown in Fig. 2. 

The wire tension was maintained by means of a 
system of pulleys and balance weights. The 
relative displacements of two symmetric points 
were assessed by the differential elongation of the 
two wires and were continuously registered. The 
different thermal expansion coefficients of the 
two different materials made possible to factorize 
out the thermal deformation of the monitoring 
system itself, subject to considerable daily and 
seasonal thermal variations under the roof 
covered by lead plates. As a result, the vault 
deformation was shown to be of about 
1 mm/yr [16]. 

The study described in [16] shows mutual 
displacements of the aforementioned points as a 
function of the monitoring time. 
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Fig. 2: Positioning of the muon telescope (T) and 
muon targets (Bl, B2, B3) on the building to be 
monitored 

It is worth remarking that the mechanical 
method actually adopted in [16] could only 
provide the measurement of the horizontal 
relative displacement of points of the wooden 
structure in opposite positions and not their 
absolute displacements with respect to a common 
reference system linked to the building structure. 

3. The proposed monitoring 
system 

3.1. Layout of the muon telescope 

Highly penetrating cosmic ray muons rain 
continuously on the Palazzo della Loggia at a rate 
of 210 000/</s, corresponding to 
10 000///(m 2 -min). The proposed monitoring 
system based on cosmic ray detection, able to 
monitor roof deformations, has the following 
general geometrical structure. Three square muon 
hodoscope modules, of about (400x400) mm 2 
area and 6.0 mm thickness, are positioned on an 
appropriate mechanical structure, axially aligned 
at a distance of 500 mm from each other. This set- 
up corresponds to the muon telescope T, shown 
Fig .3 a, having a total height of about 1 000 mm. 
The number of modules has been chosen to allow 
for a minimum level of redundancy on the 
tracking information for the crossing muon. 

The geometry of the sensitive volume of each 
hodoscope module is based on two orthogonal 
layers of 360 mm width, composed of 120 
scintillating fibers with (3 .0x3.0) mm 2 cross 
section and 400 mm length. The two layers are 
arranged along the x and y axes of a Cartesian 
reference system as shown in Fig .3b. The two 
layers provide the measurement of the crossing 
position of an incident muon in the x and y 
coordinates with a pitch of 3.0 mm. Considering a 
flat detection efficiency over the full surface of 
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the scintillating fiber, the spatial resolution on the 
coordinate of the impact is expected to be 
o = (3.0 mm)/Vl2 = 0.87 mm. 




y 



Fig. 3: (a) Structure of the muon telescope formed of 

three muon odoscope modules axially aligned at a 
distance of about 50 cm each other, and (b) sensitive 
volume of the muon odoscope module formed of two 
orthogonal layers of 120 scintillating fibers with 
(3.0x3.0) mm 2 cross section and 400 mm length. 

As shown in Fig. 2, the muon telescope T must 
be mechanically fixed to a structural element of 
the building, which is considered as the reference 
system for the measurement of the position of the 
building points to be monitored. A fourth muon 
hodoscope module, with the same geometry and 
structure as the others, must be positioned as 
muon target in each of the points to be monitored: 
in the considered case, point Bl, B2 or B3 of the 
wooden arches of the roof in Fig. 2, 
corresponding to the ones monitored by means of 
the traditional system described above. 

In this arrangement, cosmic ray muons that cross 
the full four-detectors system and the interposed 
structures of the buildings allow the continuous 
monitoring of the displacement of the muon target 
relative to the fixed muon telescope. 

The overall accuracy of the position 
measurement depends on the system geometry, 
the interposed materials, and the data collection 
time. The recognition of the signals generated by 
the same cosmic ray muon crossing the four 
detectors against background accidental events 
generated by different cosmic ray muons may be 
efficiently insured by time coincidence of the four 
detector signals [2] . 



The direction and crossing point of any cosmic 
ray traversing the muon telescope can be 
measured with an accuracy that depends on 
telescope geometry and detector granularity. The 
measured cosmic ray trajectory can be 
extrapolated from the muon telescope to the plane 
of the muon target detector, in the hypothesis that 
the trajectory of the muon is a perfect straight 
line. 

The distance between the muon telescope and 
the muon target can be deduced, with sufficient 
approximation, either from the building project 
drawings or from survey. Due to the geometry of 
the system, which limits the solid angle covered 
by cosmic ray muons traversing all four detectors, 
any uncertainty in the determination of this 
distance has negligible effects on the 
determination of the x and y coordinates of the 
extrapolated crossing point. 

3.2. Sources of measurement 
uncertainties 

The hypothesis that the muon trajectory is a 
perfect straight line is justified only in the absence 
of magnetic fields and in vacuum. The interaction 
between the terrestrial magnetic field and a 
typical 3 .0 GeV/c cosmic ray muon causes a 
deflection of about 5 /<m/m in the direction 
perpendicular to both magnetic field and particle 
velocity. Since in cosmic ray flux there is an 
excess of positive muons with respect to negative 
ones, this could induce a systematic effect on the 
measurement. However, since the measured 
displacement is obtained considering differences 
between absolute measurements performed at 
different times, the effect of the terrestrial 
magnetic field can be neglected. 

More relevant is the effect of materials 
interposed on the muon paths, which determines 
stochastic deviations of the muon trajectories due 
to multiple scattering with atomic nuclei [17, 18, 
19]. The angular deviation is proportional to the 
density of the material, to its atomic number and 
to the total amount of crossed material, and 
inversely proportional to the square of the muon 
momentum. 

The uncertainty in the prediction of the crossing 
point coordinates of the muon, in the muon target 
plane, starting from the measurement of its 
trajectory in the muon telescope, depends on 
several factors: (i) position resolution of the muon 
detectors, (ii) geometry of the muon telescope, 
(iii) distance of the muon target from the muon 
telescope, (iv) amount of multiple scattering 
angular deviations and displacement of the muon 
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trajectory, which depend on the amount of 
interposed materials and their positions. 
Being these effects mainly of stochastic nature, 
statistical distributions of the difference between 
measured crossing point coordinates in the muon 
target and the predicted crossing point coordinates 
obtained by extrapolation from the muon 
telescope is necessary to reduce the stochastic 
effects by statistical inference methods. Using the 
methods developed in [15], efficient unbiased 
estimators of the systematic displacement can be 
extracted by means of a statistical analysis of the 
distributions. The expected resolution and 
performances of the proposed measurement 
system can be analyzed by Monte Carlo 
simulations. 

Finally, it is worth noting that, depending on the 
solid angle covered by the muon telescope and on 
the position and distances of the building parts to 
be monitored, more than one muon target could 
be monitored simultaneously by the same muon 
telescope, with only some reduction of the 
acceptance of the system for the points more off- 
axis. In this way, a global and simultaneous 
stability monitoring of several parts of the 
building can be performed. 

4. Monte Carlo simulation of the 
proposed system 

The features and expected performance of the 
proposed measurement system have been 
calculated by means of a Monte Carlo simulation, 
performed by GEANT4 [20], a C++ toolkit for 
the simulation of the passage of particles through 
matter largely used in the design of nuclear and 
particle physics experiments. In the simulation, 
the geometry and the relevant structural parts of 
the Palazzo della Loggia building were taken into 
account, and the structure and component 
materials of the muon telescope and muon target 
were modeled in detail. 

In three separated simulations, the muon 
telescope was located in three different positions, 
3 .0 m below the ceiling of the large Salone 
Vanvitelliano at the first floor of the Palace, as 
shown in Fig. 2, on the vertical of each one of the 
three points of the wooden vaulted roof to be 
monitored: Bl, B2 and B3. These points are 
positioned, respectively, 0.5 m, 5.8 m, and 10.0 m 
above the ceiling of the Salone. In these three 
points a muon target was located for each of the 
three separated simulations, on the vertical of the 
corresponding muon telescope. A realistic cosmic 
ray muon generator, based on experimental data, 
was implemented in the code to simulate the 
momentum, the angular distribution, and the 



charge composition of the cosmic ray radiation at 
the sea level [21]. 

The ceiling of the Salone was modeled as a 
150 mm thick wooden layer. No vertical 
structures of the Palace were modeled, to study 
the intrinsic limits of the monitoring system 
without adding any systematic effect, nor the 
wooden vaulted roof covered by lead plates, since 
this would have introduced only a negligible 
distortion of the incoming cosmic ray spectrum 
and angular distribution, without sensitive effects 
on the system performances. 

A preliminary simulation was focused on the 
evaluation of the resolution of the muon telescope 
in the measurement of the direction of the cosmic 
ray crossing the muon telescope and in the 
prediction of the coordinates of the extrapolated 
muon crossing point on the plane of the upper 
detector module. To this aim a population of 
cosmic ray muons has been generated randomly 
on the surface of the upper detector module of the 
muon telescope. 

In each detector module the passage of the muon 
is registered when an amount of energy is 
released by ionization energy loss in one 
scintillating fiber. The measured position of the 
muon crossing point in x and y coordinates is 
defined by the position of the axis of the crossed 
scintillating fiber, on the corresponding layer, in 
the coordinate system shown in Fig. 3. 

To simplify the analysis, only the muons 
providing a single hit in all the three scintillating 
fiber layers on the same coordinate are 
considered. For this sample of cosmic ray muons, 
the three measured points on the three layers for 
each coordinate are fitted by means of a linear 
regression. The direction of the reconstructed 
straight line is compared with the direction of the 
generated muon at the entrance of the muon 
telescope and the x and y coordinates of the 
crossing point of the straight line with the upper 
surface of the upper detector module are 
compared with the coordinates of the muon origin 
point randomly extracted. 

Fig. 4a shows the distribution of the difference 
Ax between the x coordinates of a generated point 
and of the corresponding reconstructed one. Since 
the system is symmetrical, the distribution for the 
y coordinate is statistically the same as the x one. 
The distribution is obtained with a number of 
simulated events corresponding to a data taking 
time of 6 hours. The standard uncertainty in the 
determination of the muon crossing point is 
0.77 mm, compatible with the intrinsic granularity 
of the muon detector modules, which is expected 
to be 0.87 mm, as already said in section 3.1 . 
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Fig. 4: a) Distribution of the difference between the 
x coordinate of the point of generation of the cosmic 
ray on the surface of the upper detector module of 

the muon telescope and the same coordinate 
reconstructed by the straight line fit. b) Distribution 
of the difference between the projected angle on the 
x-z plane of the generated cosmic ray crossing the 
muon telescope and the same angle reconstructed 
with the straight line fit. The distributions are 
obtained with a statistic corresponding to a data 
taking time of 6 hours. 

Fig. 6b shows the distribution of the difference 
AG between the projected angle on the x-z plane of 
the generated cosmic ray and the reconstructed 
one. The standard uncertainty in the determination 
of the muon projected direction is 2.6 mrad, 
compatible with the granularity and the geometry 
of the telescope, which suggests an uncertainty of 
about 3 mrad. 

With such an accuracy in the determination of 
the direction of the cosmic ray muon, the 
expected contribution of the muon telescope 
uncertainty in the determination of the muon 
crossing point on the muon target plane 
positioned at 10 m distance is about 30 mm, that 
is, 3 mm/m. 



5. Position measurement 

uncertainty of the monitoring 
system vs data taking time 

With the proposed set-up, simulations of the 
tracking of cosmic ray muons through the 
measurement system were performed in the three 
configurations described above. For cosmic ray 
muons yielding a single hit in the four detector 
modules, the differences Ax and Ay between the 
crossing point coordinates measured in the muon 
target and the crossing point coordinates of the 
extrapolated muon trajectory measured by the 
muon telescope have been recorded. All relevant 
physical processed in the generation and tracking 
of the cosmic ray muons are taken into account in 
the simulation program. 

Cosmic ray muons generating multiple hits in 
the same detector layer are rejected from the 
analysis. Multiple hits may be produced by the 
crossing of more than one scintillating fiber in 
one layer or, alternately, by showering of the 
cosmic rays muons and possible generation of 
delta rays. The number of cosmic ray muons 
producing multiple hits in the same detector is 
about 15% of the total number entering the 
measurement system geometrical acceptance. 

Fig. 5 shows the distributions of the statistical 
variable Ax for the three configurations of the 
measurement system described above, for a data 
taking time of 15 days, corresponding to about 
31.7 10 6 cosmic ray muons crossing the muon 
target surface. The distributions of Ay are not 
shown, as they are statistically identical to the Ax 
ones. The number of events in each distribution is 
compatible with the number of expected cosmic 
ray muons entering the geometrical acceptance of 
the measurement system given by: 



■ dN A 2 A T 

d£ldA R 2 1 



Eq.l 



where d 2 NI(d£l xy dA 1 ) = 70 pi/(m 2 - s- sr) is the 
muon rate per unit area and unit solid angle 
around the zenith direction [1], Ay is the surface of 
the muon target, A 2 is the lowest surface of the 
telescope and R is their distance, and T is the 
elapsed time of the data taking in seconds. 

As the muon target and the muon telescope are 
exactly coaxial in the simulation, the Ax 
distributions are symmetric and centered at zero. 
The shape of the distributions exhibits a central 
narrow peak with very long tails on both sides. 
This shape is due both to the intrinsic uncertainty 
of the muon telescope in measuring the direction 
of the cosmic ray muon, and to the multiple 
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scattering angular deviations of the muon 
trajectories traversing the interposed materials . 

At fixed momentum and for little angles, these 
deviations follow a Gaussian law with a variance 
that depends on the inverse square of the muon 
momentum [19]. 
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Fig. 5: a) Distribution of the difference Ax between 
the crossing point coordinates detected by the muon 
target at position Bl (a), B2 (b) and B3 (c), and the 
crossing point coordinates on the muon target of the 
extrapolated muon trajectory measured by the 
muon telescope. Simulated data taking time is 15 
days. 

The latter effect dominates for larger distances 
of the muon target from the muon telescope. The 
long tails of the distributions are due in part to 
low-momentum muons, suffering larger 
deviations, and, in part, to spurious events 
corresponding to emission of delta rays or cosmic 
ray showering, most of which can be discarded 



with a more refined data analysis. At present, the 
only selection applied to these low-quality events 
is an arbitrary cut of both tails in the three 
distributions, discarding about 1.0% of the total 
events. 

RMS values of the Ax distributions are 1 .2 cm, 
6.0 cm and 8.5 cm respectively for the three 
configurations examined. They correspond to an 
estimation of the standard deviation of the parent 
populations c distr , and represent the uncertainties 
in the prediction of the crossing coordinate of the 
cosmic ray muon on the muon target detector 
extrapolated by the muon trajectory measured 
with the muon telescope. This uncertainty 
depends mainly on the distance between muon 
target and muon telescope and on the amount and 
position of the interposed materials. The mean 
value of the sample distribution represents an 
unbiased estimator of the position of the muon 
target relative to the telescope axis. The 
uncertainty on the mean value of a sample 
distribution is given by the well-known relation 
[22]: 

°mean = ° ' distr < ' E ^ 2 

where N ev is the number of events in the 
distribution. In the three configurations 
considered and for a data taking time of 15 days, 
these standard uncertainties are respectively 
0.030 mm, 0.37 mm, 0.78 mm. 

Since in the same geometrical condition N ev is 
simply proportional to the data taking time, the 
measurement standard uncertainty depends only 
on the inverse of the square root of the data taking 
time. In Fig 6 the relation of the position 
measurement standard uncertainty and the data 
taking time for the three examined conditions is 
plotted considering a data taking time up to one 
month. As time increases the measurement 
standard uncertainty decreases. By fitting the 
plots with the following general relation [15]: 

^an=C/^ Eq.3 

where C is a constant depending on the 
geometry and materials interposed and t is the 
data taking time in days, the following values for 
the constant C were obtained in the three 
conditions considered: 1.12mm- d' /2 , 

1.42mm-d' /2 , 2.96mm-d' /2 , with a standard 
uncertainty of the order of 1 % . 

As expected, the standard uncertainty of the 
measurement system depends on the geometrical 
configuration considered, since both the root 
mean square of the Ax distributions and the rate of 
useful events collected are strongly dependent on 
the geometry of the system and on the amount of 
materials interposed. Nevertheless, although 
requesting different data taking times, the 



38 



Ancona, 11-13 settembre 2014 



monitoring of the displacement of the three 
inspected points in the wooden vaulted roof of the 
Palazzo delta Loggia, using a cosmic ray tracking 
system, could provide performances compatible 
with the requested accuracy and with the time 
scale characteristic of the deformation 
phenomenon. Typical time scales, in the case of 
Palazzo della Loggia and, in general, for 
historical buildings, may span over several years. 




time (days) 



Fig. 6: Standard measurement uncertainty as a 
function of data taking time 

In position Bl, a measurement standard 
uncertainty of the order of 0.1 mm may be 
achieved in about one day of data taking, whereas 
a standard uncertainty of the order of 0.5 mm may 
be achieved in a week of data taking in position 
B2 and in one month of data taking in position 
B3, where muon target and muon telescope are 
positioned 13.0 m far apart. 

6. Performances of the proposed 
system vs the traditional 
approach 

To demonstrate how a monitoring system based 
on cosmic ray muon tracking may perform better 
than or well as more traditional systems, provided 
that the acquisition time is long enough, the 
results shown in Fig. 6 have been compared to 
those performed with the traditional monitoring 
system described in section 2.2 and presented in 
[16]. In [16] both the cyclic seasonal deformation 
of few millimeters for a position of the wooden 
vaulted roof structure of the Palazzo della Loggia 
and the general trend of displacement with a rate 
of deformation of the order of 1-mm/yr are 
shown. 

A first result is that a monitoring system based 
on cosmic ray tracking would be largely sensitive 
enough to detect the progressive deflection of 1 
mm/yr of the vaulted wooden roof of the Palazzo 
della Loggia. Indeed, looking at Fig. 6, it is 
possible to deduce that, in the worst studied case, 
the proposed simulated system reaches a standard 
measurement uncertainty of 1 mm after an 
elapsed time of about 10 days. This means that 



the system is able to identify a displacement of 
1 mm, with a 68.3% confidence interval, in 
10 days. 

Moreover, in [16], the measured seasonal 
displacement of the roof structure on a period of 
several years are reported. To investigate how a 
monitoring system based on cosmic ray muon 
tracking could follow cyclic seasonal 
deformations as well as systematic ones, a 
seasonal displacement of the three positions Bl, 
B2 and B3 in Fig. 2 was reproduced in the 
simulation program as a function of time. The 
behavior of the structure in point B2 for the first 
year of data taking, experimentally measured 
thanks to the traditional mechanical monitoring 
system previously described, was adopted as a 
realistic model of the phenomenon. A cosmic ray 
data taking one year long was simulated to 
reproduce the measurement of the possible 
seasonal displacements of the roof structure in 
points Bl, B2 and B3 by the proposed monitoring 
system. 

In Fig. 7, the displacement with time 
corresponding to the supposed seasonal 
deformation (the same for the three points) is 
shown as a continuous line. The results of the 
simulated measurements using the sample mean 
of the position of the muon target, displaced 
following the assumed structural deformation, is 
shown, with sampling rate of one week, two 
weeks and one month respectively for points B 1 , 
B2andB3. 

The ability of the proposed measurement system 
to follow seasonal displacements of few 
millimeters and, even more so, also systematic 
ones is demonstrated. 

7. Discussion and conclusion 

The comparison between the numerical results 
of the proposed feasibility study, presented in 
Fig. 6, and the experimental results reported in 
[16], demonstrates that a measurement system 
based on cosmic ray muon tracking is suitable to 
monitor the static stability of historical buildings. 
In the particular case study considered, seasonal 
displacements can be identified, and systematic 
displacements that could affect building stability 
can be measured with a standard uncertainty equal 
to or even better than the one associated with 
traditional methods. 

The assessment of the overall monitoring system 
standard uncertainty could be significantly 
improved both by modifying some geometrical 
parameters of the system and improving the data 
analysis: 
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Fig. 7: a) Seasonal deformation of few millimeters 
assumed for the point at position Bl Fig. 2. Dots 
indicate the simulated measurement of the position 
of the muon target with the muon stability 
monitoring system using the sample mean value, 
with sampling rate of one week. The measurement 
uncertainty is given by equation Eq.2. b) The same 
for the point in position B2, with sampling rate of 
two weeks, c) The same for point in position B3, 
with sampling rate of one month. 

1. In [15] it was demonstrated that a more 
efficient unbiased estimator of the position 
and possible displacement of the muon 
target detector can be obtained by fitting the 
parent population of the statistical variables 
Ax and Ay by an appropriate analytical 
function, and by determining their shape 
parameters using a best fit procedure. In this 
way, the standard uncertainty is expected to 
be improved by a factor of 2 to 3 [13] . 



2. The same work [15] points out that the low 
momentum muons can be discarded by the 
data sample by means of a muon absorber in 
the form of a thick iron plate positioned 
below the muon telescope and followed by a 
further plane scintillation counter in 
coincidence. By eliminating the low 
momentum muons, which suffer the largest 
multiple scattering deviations, and the 
possible electron showering that adds 
spurious events to the sample, the long tails 
of the Ax and Ay distribution can be 
reduced, thus reducing the root mean square 
Odistr of the distributions. 

3. A possible way to improve the accuracy of 
the monitoring system and to decrease data 
collection time involves both system 
parameters and data analysis. Instead of 
using only one detector module, the muon 
target could consist of a couple of particle 
detector modules assembled to form a 
telescope, allowing both the crossing 
position and the direction of the incident 
muon cosmic ray to be measured at the 
muon target point. 

In this way, when the multiple scattering 
deviation is concentrated in only one main 
layer (an interposed floor for example), the 
deviation effect on the muon trajectories can 
be largely compensated and it was estimated 
by Monte Carlo simulation that the system 
measurement uncertainty previously reported 
may be reduced up to a factor 3 . 
When the building structures interposed 
between the muon telescope and the muon 
target are more complex (presence of more 
than one floor or vertical walls), a 
reconstruction algorithm based on a 
likelihood function and on Monte Carlo 
simulation of the muon statistical behavior, 
can be implemented to determine the more 
probable position of the muon target 
accounting for the known distribution of the 
interposed structures and composing 
materials . 

4. Finally, an increment of 50% of the single 
module size corresponds to a reduction of a 
factor 2.2 of the data collection time needed 
to obtain the same standard uncertainty on 
the position measurement. This solution has 
to be compatible principally with the cost of 
suitable particle detectors and with the 
complexity of their control system. 

In conclusion, cosmic ray muon detection 
techniques have been investigated for 
measurement applications in the field of civil 
engineering and have been demonstrated to be 
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particularly suitable for static monitoring of 
historical buildings, where the evolution of the 
deformation phenomena is of the order of months 
or years. Appealing features of the proposed 
monitoring system are: (i) the use of a natural and 
ubiquitous source of radiation avoiding any 
problem of radiation protection; (ii) its 
applicability also in presence of horizontal and/or 
vertical building structures interposed between the 
reference system and the parts to be monitored; 
(iii) the limited invasiveness, and the flexibility 
and ease of installation of the monitoring system 
devices; (iv) the possibility to design a global 
monitoring system, where the position of different 
points of the building may be simultaneously 
monitored relative to the same reference system; 
(v) the use of well-known physical principles and 
established technologies in the field of nuclear 
and particle physics to build up the cosmic ray 
muon detectors featuring the characteristics 
suitable to satisfy the requirements of any specific 
application. 
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Abstract 

This paper aims at describing the dynamic 
behaviour of polymeric shape-memory stents. In 
particular, considering the recovering phase, 
measurements have been performed to find a 
relation between different configurations and 
forces, which the specimens exerts. Shape memory 
materials, in fact, can exert a force when 
recovering their permanent form. 
An analytic model that describes the shape of an 
expanding specimen, based on a vision system, 
has been developed and experimentally 
calibrated. 

A measurement device has been specially 
designed to measure forces, which the expanding 
specimen exerts . This device is based on six 
trolleys free to fall, frictionless, thanks to 
pneumostatic bearing, on an inclined rail. So 
entity of the force applied on a shape memory 
specimen, constant along the whole movement, is 
determined by the inclination of the guide. 

1. Introduction 

With the name "shape-memory polymers" 
(SMPs) literature refers to a category of smart 
materials, which have the ability to be deformed 
from a "permanent" shape to a "temporary" one, 
and to revert back to the "permanent" shape once 
an external stimulus, usually a change in 
temperature, is applied. 

Thanks to the peculiar characteristics of 
polymers, SMPs may offer various advantages 
against their metallic counterparts (shape memory 
alloys) [2, 3], such as being lightweight and easy 
to process, being capable of recovering large 
strains and a relatively easily tailoring of the 
temperature that triggers the shape -memory 
effect. Literature is investigating the possibilities 
to employ SMPs in realistic situations, and in 
particular in the fields of sensors, actuators, 
biomaterials and smart textiles [2, 4-11]. It is, 
however, to be remarked that shape memory 
polymers are typically not capable of exerting 
high levels of stress, so that the capabilities of 
these systems meet the requirements of what can 



be considered a soft actuation, which could be 
employed for biomedical applications, due to the 
compliance of soft tissues. 

This paper investigates the capability of 
polymer-based star-shaped specimens to behave 
as actuators. In particular, the relation between 
their shape changes during the recovering phase 
and the forces that they can exert. 

The activity of measurement described in this 
paper regards two different mechanical quantities: 
force and deformation/opening level. 

For the measure of force a specific device has 
been designed and realized, while for the 
assessment of the deformation/opening level an 
image -based method has been proposed. 

2. Experimental kinematics 

The objects of the experimental investigations 
were actuators of a semicrystalline crosslinked 
polymer, poly (e-caprolactone) (PCL), which 
were prepared with a tubular shape as the 
"permanent" one. The choice of PCL networks as 
shape memory systems for our analysis is based 
on the fact that this class of polymers may be of 
large interest for biomedical applications and 
since their transformation temperature can be 
finely tuned by means of a proper choice of the 
material crossllink density, as shown in previous 
papers [16, 17]. The crosslinked structure was 
obtained starting from hydroxyl-terminated 
poly(s-caprolactone) (PCL), with a molecular 
weight of 7000 g mol-1 , later suitably modified to 
obtain methacrylate -terminated PCL precursors 
and crosslinked by a free-radical thermally 
activated process in presence of dicumyl peroxide 
as radical initiator (at 105°C for 5 hours). 

The tubular specimens were prepared by pouring 
the reactive mixture in silicon molds with a cavity 
having the shape of a tube (height = 40 mm, outer 
diameter = 25 mm, thickness =1.5 mm). Special 
attention was taken to reduce the presence of air 
bubbles that may be trapped in the cross-linked 
materials, maintaining the melted methacrylate- 
terminated polymers under reduced pressure for 
several hours before the thermal curing; more 
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details concerning the preparation of these 
materials are reported in a companion paper [18]. 

Differential scanning calorimetry (DSC) tests 
carried out on the material obtained revealed the 
presence of a semicrystalline structure 
(crystallinity content: about 40%), with a melting 
temperature Tm = 53°C and a crystallization 
temperature Tc = 18°C. 

The shape memory response of the material was 
tested by first setting the specimen to a 
"temporary" configuration. This step, typically 
called "programming", was carried out by heating 
the cylinder at a about 80 °C (i.e. above Tm) and, 
by employing an ad-hoc designed crimping tool 
(Figure 1), which allowed to fold them in a six 
arms star shaped configuration. 




Fig. 1: i) side-view of the crimping fixture; ii) top- 
view of the crimping fixture, without the top closure 
plate, and picture of an undeformed (A) and 
partially crimped (B) specimen. 

Once crimped, the specimen was cooled at about 
-17 °C (i.e. well below Tc) for at least 10 minutes 
while maintaining the crimping arms blocked, so 
to fix this compact configuration. The folded 
specimens were than subjected to isothermal 
stress-free recovery test by immersion into a 
water bath, previously heated at various 
temperatures, chosen from room temperature 
(20 °C) to above the material melting temperature 
(53 °C). The recovery process was monitored as a 
function of time by means of a photocamera with 
a sampling frequency of 1 frame/s. In this case a 
color image (32 bit 4288x2848 pixel) is taken 
from by a Nikon D70 camera. 

The star shaped actuator is monitored during a 
period of 17 seconds. A photograph of the 
opening stent is realized at the starting time and, 
as already said, every second during a 
transformation period. During this lapse of time 
the specimen receives thermal energy and realizes 
a conversion to mechanical energy through a 
shape change, as shown in Figure 2. In this 
experimental test, the heat transfer to the actuator 
is performed with a liquid media, where the 
actuator is completely immersed. 

To describe the kinematics of the stent during its 
motion, an analytic model has been assumed and 
experimentally calibrated. 



Fig. 2: Transition from temporary shape to 
permanent shape 

In the initial instant, the shape could be 
described as a linear combination of rhodonea 
functions, and in the final instant as an annulus. 
The proposed model linearly combines these 
functions with parameters that are time- 
dependent. Figure 3 shows how the model (blue 
line) fits experimental data (red line), determined 
using an edge detection and filtering algorithm, 
specially developed. 



Fig. 3: Comparison between the detected boundary 
(red line) and the model boundary (blue line) 

A complete description of the mathematical 
approach has been described in a dedicated paper 
presented by the authors [15]. 

The same approach has been used for the 
identification of the material opening 
level/deformation during the force controlled test, 
performed in order to find a relation between 
shape at a given opening level and the 
corresponding exerted force. During the force- 
deformation test the fluid media used for heating 
the polymer-based star-shaped specimens was the 
air at a controlled temperature. 

For the opening level an evaluation of the 
uncertainty has been performed considering the 
results of the calibration of the kinematics 
mathematical model described in [15]. 

Analyzing the results of the model calibration 
and considering the resolution and the 
characteristic of the vision system used (HR 
digital camera) a value for the standard 
uncertainty, considered in further tests as "B" 
type, is 0.05 mm. 

3. Force measuring device 

A measuring device, shown in Figure 4, has 
been specifically developed to measure the 
radial force [12-14] that each specimen exerts 
during the transition from temporary to permanent 
shape. This transition takes place in air at a 
uniform and controlled trigger temperature. The 
device consists of six trolleys, in radial position 
with respect to the specimen. Each one is free to 
fall, frictionless thanks to pneumostatic bearing, 
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on an inclined rail. The angle (0) of the rail can be 
modified, to vary the load (L) applied to the 
specimen. The measuring device is based on null 
measurement method and the radial force (F stent ) 
can be unequivocally determined when the 
component of the load of the trolley results equal 
to the component of the force which the opening 
specimen exerts, both components directed along 
the inclined rail. This is a static situation and the 
specimen keeps its shape unchanged. 




Fig. 4: measuring device 



The simplest model for the evaluation of the 
F stent can be express by the relation: 

where m represents the value of the mass of a 
single trolley, g is the gravity acceleration, 6 is the 
angle of the rail and F friction represents the residual 
part of friction between rail and trolley. 

In the model of the measurement the friction 
force has been modelled as a statistical quantity 
with null average ad standard deviation of 
0.010 N. More investigation in that regard is 
needed and will be performed in later stages of 
this project using a comparison dynamometer. 

From Eq. 1, a first evaluation of the 
measurement uncertainty (cat. B) can be 
calculated using the ISO 13005 technique as 
shown in tables: 
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with: 

^ stem = 0.275 N and «(0 = 0.014 N 
therefore the relative uncertainty is assessed at 
5.2% and the main uncertainty contribution is the 
presence of friction between rail and trolley. 

Similar results can also be obtained with other 
rail's inclination. 



4. Mechanical characterization 

When the equilibrium between the 
aforementioned components is achieved, F stent can 
be assessed using equations of static equilibrium 
of the system. Figure 5 shows some static 
configurations and the corresponding force 
exerted by the specimen. 



0=7° 0=5° 0=2° 




Fig. 5: Each configuration of the expanding 
specimen corresponds to a radial force exerted by 
the specimen itself 

The purpose of the experimental activity is to 
characterize mechanical behaviour of shape- 
memory materials, in particular to determine a 
possible relation between forces exerted by a 
specimen made of the considered material and its 
geometrical shape during the recovering phase, 
when the specimen is expanding. Forces are 
determined for different external load imposed 
and temperature conditions. In this way, it is 
possible to verify the actual usability of this kind 
of materials in biomedical applications . 

The experimental activity allowed to 
characterize the main behaviour of the material. 
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The first indication is that the material has a 
variation of response in different examination as it 
is shown in next figure 6. 

Figure 6 shows the level of opening of the 
specimen with a constant force level applied (rail- 
trolley inclination) in different testing. It is 
possible to observe that the first testing allows to 
the materials a greater level of expansion than the 
others tests. Each test in Figure 6 is performed in 
the same temperature and load conditions, using 
the same specimen. The variability between 
different test for the same specimen in the same 
condition testing (temperature, external load) is 
greater then uncertainty. 

A different test has been also performed with a 
set of different specimens of the same material at 
their first test. For each experiment a different 
level of external load (rail-trolley inclination) has 
been applied. 




Fig. 6: repeatability 

Results obtained show how the shape memory 
material works: when it is compressed to the 
non permanent shape (small entity of opening 
level) the force that the specimen can exert is the 
highest. During the expansion, this force 
decreases quickly and it became null at the 
maximum expansion. 



heights of specimens, in order to verify how the 
mechanical behavior of the devices is influenced 
by these variables . 

5. Developments 

Tests carried out have shown how the device 
developed for the force measurement allows to 
test only specimens of shape memory material 
capable of exercising limited forces up to 0.55 N. 
In order to increase the full scale of the force 
measuring instrument it will be necessary to 
install calibrated weights on the singular trolley 
for increase its mass. 

A second activity, which must be further 
developed, is the direct metrological properties 
verification of the dynamometer proposed. In 
particular, the residual friction force entity should 
be assessed, also in relation to the increase in 
weight of the trolleys . 

6. Conclusions 

The present work concerns a methodology for 
the mechanical characterization of shape memory 
materials. 

The proposed procedure requires the use of a 
vision system for assessing the level of shape 
changes (opening level). Forces exerted during 
the expansion of tested specimens were evaluated 
with a measuring tool specifically designed and 
developed. 

The difficulty of the force measurement consists 
in the application of a load of very low intensity 
independently from its point of application. The 
solution adopted is the use of a trolley free to fall 
along an inclined rail with a pneumatic bearing. 

Tests carried out have allowed to define a force- 
deformation curve that allows to evaluate material 
performances for its use as a stent in biomedical 
applications. 
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Fig. 7: mechanical behavior of the material 

Many tests were conducted with various material 
compositions, different thicknesses and different 
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Abstract 

In planetary exploration space missions, motion 
measurement of a vehicle on the surface of a 
planet is a very important task. The main 
contribution of this work is an experimental 
comparison among visual odometry systems using 
lenses with three different focal lengths (an ultra 
wide angle, a medium wide angle and a telephoto 
lens). For each focal length a complete 
calibration procedure is performed using a 
suitable reference instrument and taking into 
account different types of motion. The calibration 
procedure comprises a rigorous uncertainty 
analysis. Experimentally obtained uncertainties 
are compared, pointing out strengths and 
weaknesses of employing wide angle or telephoto 
lenses for motion measurement. 

1. Introduction 

In planetary exploration, motion measurement of 
a vehicle on the surface of a planet should be very 
accurate in order to track the vehicle also for long 
paths. The odometric evaluation of vehicle 
position and attitude performed measuring the 
rotation of wheels has wide uncertainty due to 
slippage of wheels on a natural, often sandy or 
slippery, surface. Moreover, on extraterrestrial 
planets GPS -like positioning systems are not yet 
available and inertial navigation sensors exhibit 
unacceptable drifts. Thus, the need of a reliable 
and accurate motion instrument is particularly 
relevant. In this work a vision-based instrument 
for displacement and rotation measurement is 
calibrated using a high precision motor-driven 
rotary stage and a linear slide as reference 
instruments. To evaluate the measurement 
performance along paths longer than the reference 
linear slide, the visual odometry system is also 
mounted on a laboratory vehicle, which is driven 
along a closed path and the whole travelled 
trajectory is acquired by the stereo vision system. 

The use of stereo systems for visual-odometry 
(VO) is well known, e.g. see [1]-[12], but is a 
still-open research subject, as proved by many 
recent papers, e.g. see [13] -[16]. An interesting 



overview and an introducing tutorial on visual 
odometry can be found in [11] and [12]. In visual 
odometry, the displacement and rotation of a 
stereo vision system are measured through the 
images taken by the two cameras. Stereo- 
processing allows estimation of the three 
dimensional (3D) location of landmarks observed 
by a stereo-camera. If the same landmarks are 
acquired and detected by a moving stereo-system 
in two subsequent positions, the two 3D point 
clouds allow to evaluate the stereo camera 
movement (position and orientation) that took 
place between the acquisitions of the two stereo 
images. The whole trajectory is then calculated 
combining each motion step. 

Particularly, reference [5] describes a method for 
visual odometry based on a stereo camera, 
providing experimental results gathered during 
development and flight phases of the NASA's 
twin Mars Exploration Rovers Spirit and 
Opportunity landed on the surface of Mars in 
January 2004. Reference [6], for the first time, 
analyzes the need of an anisotropic uncertainty 
modeling for 3D acquired landmarks. Reference 
[9] emphasizes the importance of a detailed and 
correct uncertainty evaluation, and describes a 
method for visual-odometry that allows to reduce 
the final measurement uncertainty. The methods 
described in [1] - [10] require that the same 3D 
landmarks are observed in two subsequent 
acquired stereo images; this requirement limits 
the maximum translation and, even more, the 
maximum rotation that can be measured in one 
step. To follow a whole trajectory, e.g. of a 
vehicle with the stereo system mounted, several 
incremental measurements are performed and 
combined together, with the warning that 
uncertainty may significantly increase with the 
number of steps, i.e. small errors in each 
measurement step can eventually cause large 
errors in the estimated trajectory. This drawback 
can be solved using more global approaches (e.g. 
Simultaneous Localization And Mapping, loop- 
closing, semi-global optimization). However, in 
the present work the attention is focused on 
comparing the measurement uncertainties 
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obtained by a visual odometry system using 
lenses with different focal lengths. The main 
purpose is to experimentally compare the 
metrological behavior obtained with different 
lenses, and not to find the best position 
measurement method based on a vision system. 
Thus, more global approaches are beyond the 
purpose of this work and are not taken into 
account. 

The main aim of the presented comparison 
among different focal lengths and, thus, different 
values of field of view (FOV) is to yield useful 
advices to the problem of lens selection for the 
stereo camera of a visual odometry system. 
Reference [4] performed numerical simulations to 
evaluate which FOV yields the smallest long- 
range error of a visual odometry system. Authors 
analyzed several values of FOV between 15° and 
90° using a numerical simulation procedure, and 
found that the optimal FOV is approximately 35°. 
In [4], the obtained numerical results are not 
validated by an experimental comparison. Despite 
the optimal values of FOV found in [4] , reference 
[5] says that their visual odometry algorithm has 
been experimentally tested using the Jet 
Propulsion Laboratory's Rocky 8 rover, whose 
cameras have a horizontal FOV equal to 80° and 
64° along vertical direction. Other tests were run 
on the MER Surface System Testbed Lite rover 
with a 120° FOV cameras and, of course, on the 
Mars rovers with a 45° FOV. However, in [5] 
there is not a rigorous comparison among 
experimental results obtained with these different 
FOV values. Moreover, different values of FOV 
are associated with different cameras and 
different rovers, thus, a direct comparison could 
not be meaningful, since FOV is not the only 
parameter that changes from one case to another 
one. Our contribution is a direct experimental 
comparison among visual odometry systems with 
three different focal lengths and FOV (an ultra 
wide angle, a medium wide angle and a telephoto 
lens), performed with all other influencing 
parameters kept constant (same cameras, same 
relative positions of cameras, same elevation 
angle of cameras, same imposed rotary and linear 
motions) and with a rigorous uncertainty analysis 
according to [17]-[18]. 

In [28] , we applied about the same VO algorithm 
described in the present manuscript to the 
measurement of a vehicle trajectory. The vehicle 
was driven along a closed path and was brought to 
the initial position with a positioning uncertainty 
of about 1cm. The on-board mounted stereo 
camera acquired two video sequences with fixed 
frame rate during the motion. The optical 
encoders mounted on the wheels were not able to 



acquire the trajectory with an uncertainty one 
order of magnitude better than the VO system. 
Thus, only the final position could be used to 
analyze the errors and the uncertainty achieved by 
the VO system. With reference to [28], the major 
contribution of the present manuscript is to 
describe the behavior of different VO systems 
using a complete new laboratory set-up, which 
allows to analyze the errors and uncertainties 
obtained by each VO system for all the acquired 
motion steps and not only at the end of the 
motion. With the new experimental set-up, we 
have the possibility to measure both linear 
displacements and rotations imposed to the stereo 
camera during all the motion steps with an 
uncertainty an order of magnitude better than that 
of the VO system. Thus a more detailed analysis 
can be performed than in [28] . 

In section 2, the paper describes a measurement 
method broadly based on the NASA rovers 
approach [4], [5] and slightly updated combining 
together subroutines and procedures more recent 
and advanced. Section 3 presents the performed 
uncertainty analysis, and section 4 discusses the 
experimental set-up and the obtained results . 

2. Measurement algorithm 

In this section, the procedure employed to 
perform the displacement and rotation 
measurement of the stereo system is described. In 
each motion step, the goal is to calculate the 
displacement and rotation of a calibrated stereo 
system using the images acquired in two 
subsequent positions. The whole trajectory of the 
vehicle is then evaluated combining each single 
motion step. Rotation is described by a sequence 
of Euler angles around axes X, Y, Z. 

The procedure begins with the detection of 
image features (keypoints) which are the 
projections of physical landmarks in the two 
cameras. After the intrinsic and extrinsic 
parameters of a stereo system are carefully 
determined (the stereo camera is calibrated as 
described in [27]), one of the main uncertainty 
sources is the position on the image plane of the 
features detected and matched in corresponding 
images. Several feature detectors and descriptors 
were considered in the preliminary phases of this 
work, see [19]-[26]. Particularly, reference [23] 
compares several different detectors invariant to 
scale and affine transformations and finds out that 
the best results are obtained by the Hessian- Affine 
detector and the Maximally Stable Extremal 
Regions (MSER) approach. These methods are a 
good choice to find features even in presence of 
relatively wide scale and/or rotation variations 
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and of affine transformations between images. 
However, in our work we used two well known 
algorithms: the Harris detector [19] and the Scale 
Invariant Feature Transform (SIFT) detector [22] , 
since two selected approaches are widely used 
and we do not need to cope with wide affine 
transformations when we apply the VO algorithm 
to subsequent images. 

Reference [26] makes a comparison among 
different descriptors and concludes that the SIFT 
and the Gradient Location and Orientation 
Histogram GLOH descriptors are the best. For 
this reason, the SIFT descriptor is chosen to 
describe the image neighborhood of each detected 
feature and to reliably perform feature matching 
between stereo pairs and between corresponding 
features in subsequent vehicle movements. The 
SIFT descriptor is used for all features: those 
detected by SIFT detector and also those found by 
the Harris detector. This choice is one of the 
major difference between the employed algorithm 
and that of NASA rovers [4] -[5], which employ a 
Forstner or Harris operator as feature detector, 
and rely upon cross correlation along epipolar 
lines for feature matching. This is also a 
difference with our previous work [28]: in the 
present application we found out that the Hessian- 
Affine detector does not exhibit a significant 
advantage over the combination of Harris and 
SIFT detectors. 

After the corresponding 2D features are 
detected, a triangulation phase allows to compute 
the 3D coordinates of the physical landmarks 
acquired by the cameras. The middle point 
algorithm is used for triangulation, as in [5], for 
more details see [10]. 

The uncertainties of 3D points (their 3x3 
covariance matrix) are needed for the evaluation 
of stereo system displacement and rotation as 
described in subsection 2.1. Uncertainty 
evaluation for triangulated 3D points becomes an 
uncertainty propagation task, which is performed 
by the Kline-McClintock formula, see GUM [17]. 
This method is selected, instead of the Monte 
Carlo propagation approach, since the calculation 
is embedded in the algorithm that calculate the 
displacement and rotation of the stereo camera, 
and is performed for all detected features (see 
subsection 2.1). Thus, the use of a Monte Carlo 
simulation could lead to unacceptable time delays 
in the position and attitude evaluation of the 
stereo camera. 



2.1. Stereo system displacement and 
rotation 

When the stereo system is calibrated (intrinsic 
parameters of both cameras and the position and 
orientation of camera 2 with reference to camera 
1 are known), the 3D points Pl X can be calculated 
for all features detected by both cameras when the 
vision system is in an initial position PI. The 
notation PI X means that these points are expressed 
in the reference frame attached to the first camera, 
when the vision system is in the initial position 
PI. When the vision system is moved (cameras 
are rigidly connected) from the initial position PI 
to a subsequent position P2, the same procedure 
can be used to compute the 3D vectors ra X of the 
same features detected by both cameras in the 
second position P2 and expressed in the new 
frame 1 attached to the first camera. For each 
feature that is detected by both cameras in both 
positions PI, P2, the following equation can be 
written: 



Eq.l 



Where pi** is the rotation matrix from frame 1 
in the second position P2 to frame 1 in the initial 
position PI; p1 Pp2,pi is the origin of frame 1 in P2 
with reference to the origin of frame 1 in PI and 
expressed in PI; i= 1, n, with n being the 
number of common detected and matched 
features. Vector Pl ~P P2? pi and the Euler angles that 

define matrix />2 R are the numerical output 
values of the whole measurement procedure and, 
in this work are evaluated in two steps in a similar 
way as in [5]: first, a less accurate motion is 
estimated by least squares estimation embedded 
within a random sample consensus (RANSAC) 
process to remove outliers; then, a maximum 
likelihood motion estimation is performed 
minimizing a non linear problem. The main 
differences with [5] are: the feature detection and 
matching algorithms employed in this work are 
more advanced; the linear least squares approach 
is derived from [7]; the non linear minimization 
procedure comprises the Levenberg-Marquardt 
algorithm. 

In the linear motion evaluation an error vector e, 
is defined for each couple of 3D points P1 X, P2 X 
detected in both the first PI and the second P2 
position and then the cost function E to be 
minimized is calculated: 



l P2,P\ 



Eq.2 



50 



Ancona, 11-13 settembre 2014 



E = ihf E <i-3 

i=i 

In order to separate the evaluation of rotation 
and translation, the centers of the two point clouds 
are subtracted from the 3D points. Using the 
formulas explained in [7] , first the rotation matrix 
is evaluated and then the translation pi Pp2,pi of the 
stereo system is calculated rearranging Eq.l. 

In the non linear step also the covariance 
matrices of the 3D landmarks are taken into 
account. For the non linear analysis, the errors 
defined for each 3D feature in Eq.2 is weighted 
using the inverse of a combined covariance 
matrix. For each 3D point evaluated with the 
stereo system in the first position, the 
corresponding covariance matrix PI Q is 
evaluated; then, this evaluation is repeated for the 
same feature in the second position to obtain the 
covariance matrix P2 Q. P1 C, and P2 C, of the same 
3D landmark are combined together and the new 
cost function E nl (minimized by the Levenberg- 
Marquardt algorithm) is: 

C,. = «C j + P P 2R- P2 Cr^R r Eq.4 
Em=Yj?i -Cr'-e,-) Eq.5 

i=\ 

In this way each component of the error vector e,- 
of the feature i can have a different weight that 
takes into account the uncertainty of feature i 
along the considered direction. This 
heteroscedastic modeling of uncertainty, i.e. 
inhomogeneous (it may be different from point to 
point) and anisotropic (it may be different in each 
direction) is particularly useful in stereo systems 
whose baseline (i.e. distance between cameras) is 
small and the uncertainty along a direction 
roughly parallel to the optical axes may be much 
greater than in the other directions. A properly 
tuned RANSAC algorithm used in the linear 
phase of motion estimation allows to identify and 
to exclude possible outlier landmarks. Thus, the 
non linear phase is performed taking into account 
only 3D points that passed RANSAC algorithm. 

3. Uncertainty analysis 

A detailed uncertainty analysis is carried out 
according to the metrological procedures 
described in [17] and [18]. The position and 
orientation uncertainty of the stereo camera is 
evaluated as in an indirect measurement by a 
Monte Carlo propagation approach. Along the 
imposed motion, the position and orientation of 
the stereo camera are obtained combining every 



motion step and are treated as the output 
quantities of an indirect measurement. This 
uncertainty analysis can be performed off-line 
after the trajectory has been evaluated, thus, the 
aim is to evaluate uncertainties of position and 
attitude of the stereo system as precisely as 
possible, and not to calculate them in a fast way. 
Moreover, the whole measurement algorithm is 
highly non linear. Thus, a Monte Carlo simulation 
is employed for this off-line uncertainty analysis 
instead of the propagation formula embedded in 
the calculation algorithm as described in section 2 
and that is an approximated method. In the new 
set-up, all the translations and/or rotations 
imposed to the stereo camera are known with an 
uncertainty at least one order of magnitude better 
than the measurements obtained by the VO 
system. Thus, the compatibility of the 
experimentally obtained errors (the difference 
between the values obtained by the VO system 
and those obtained by the reference 
instrumentation) with the uncertainties evaluated 
by the Monte Carlo method can be checked for all 
the positions and orientations imposed. 

Several uncertainty sources are analyzed and 
evaluated using experimental tests. Particularly, 
the uncertainty associated with the following 
quantities are taken into account: intrinsic and 
extrinsic parameters of the stereo system, whose 
uncertainties are evaluated using the camera 
calibration procedure described in [27], and the 
positions of image features, whose uncertainty is 
affected by several contributions. In this paper 
two main contributions are evaluated for image 
features: the reading uncertainty (image noise) of 
the sensors, and the residual uncorrected optical 
distortions that the distortion model is not able to 
remove. According to [17] and [18], all 
uncertainty sources are expressed by a probability 
density function (PDF) and are then propagated to 
the output displacement of the stereo system using 
a Monte Carlo simulation. In this way, for each 
motion step the corresponding position and 
orientation uncertainties. 

4. Experimetal set-up and Results 

In the calibration procedure of the visual 
odometry system, the stereo camera is mounted 
both on a high precision motor-driven rotary stage 
and on a linear slide. The rotary stage is driven by 
a stepper motor capable of a resolution equal to 
1.125 arcseconds, while the linear slide is 
provided with a graduated scale, in order to 
compare the measurements obtained by the visual 
system respectively with known rotation angles 
and with known linear displacements. Fig.l 
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depicts the employed set-up comprising the stereo 
camera mounted on the rotary stage, which can 
translate along the linear slide. 




Fig. 1: Experimental set-up 



In [4], the simulated stereo camera had a fixed 
baseline equal to 10cm, a height from ground of 
1.4 m and a downward tilt of 30°. The rover 
movement was simulated along a 500m path and 
with simulated measurement steps ranging 
between 30 and 70 cm. Reference [4] says that the 
optimal FOV remained between 30° and 40° also 
changing the length of measurement steps. In our 
laboratory set-up, we use a stereo camera having a 
baseline of 54 cm, a height from ground of 1.1 m 
and a downward tilt angle of 0° . Three lenses with 
different focal lengths are employed: 6mm, 10.4 
mm, 50 mm. These nominal lengths, combined 
with the sensor size of the cameras, yield 
respectively the following FOV values: 86° x 53°, 
57° x 32°, 13° x 7°. According to [4], we expect 
that the medium wide angle lens (57° x 32° FOV) 
should be the best choice, while the ultra wide 
angle lens (86° x 53° FOV), that exhibits the 
largest distance from the numerically evaluated 
optimal range (30° -40°), should be worst case. 

In our lab, a translation test with a total travel of 
1350mm using the linear slide and a rotation test 
from 0° to 90° using the rotary stage are 
performed. In the translation test, the stereo 
camera is incrementally translated along the linear 
slide with a motion step of 50 mm. After each 
step the motion is stopped and two images are 
acquired. Thus, two stereo images are available 
exactly every 50 mm. The same procedure is 
repeated in the rotation test using the rotary stage. 
In this case the stereo images are acquired every 
1°. 

The results obtained using the three lenses are 
summarized by the following figures. Fig. 2- Fig. 
5 show results obtained in the linear translation 
test with the linear slide, while Fig. 6 - Fig. 9 deal 



with the rotation test using the rotational stage. 
Fig. 2 shows the total displacement measured by 
the VO system vs. the one imposed along the 
linear slide. Fig. 3 illustrates the measured total 
errors for each motion step, which are the 
differences between the measured and the 
imposed total displacements. 
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Fig. 2: Measured total displacement in the linear 
translation test 
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Fig. 3: Total displacement error in the linear 

translation test 
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Fig. 4: Total position standard uncertainty along x 
axis in the linear translation test 
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Position standard uncertainty along z axis 
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Fig. 5: Total position standard uncertainty along z 
axis in the linear translation test 



Fig. 4 and Fig. 5 depict the standard 
uncertainties evaluated by the Monte Carlo 
approach for the total position along respectively 
the x and z axis (z axis is parallel to the linear 
slide, while the x axis is orthogonal). These 
uncertainties are obtained combining together all 
the single uncertainties of each motion step from 
the first one to the considered position. 

Results slightly different from the expected ones 
emerge from the experimental tests. From Fig. 2 - 
Fig. 5, it is clear that all three lenses are able to 
measure the linear translation along a direction 
parallel to the optical axes of the stereo cameras. 
However, the worst case is the telephoto lens 
(f=50mm) and not the ultra wide angle one 
(f=6mm), as expected from [4]. In this linear test, 
the telephoto lens is always the worst case among 
the three lenses, but it behaves in a better way 
than during the tests described in [28]. 

In the translation test, both the two wide angle 
lenses perform very well, with similar results. In 
Fig. 3 the measurement errors obtained by 
medium wide angle (f=10mm) are lower than 
those obtained by the ultra wide angle except than 
in the first part of the translation (<400 mm). 
Looking at the evaluated uncertainties, the ultra 
wide angle lens achieves better values along the z 
axis than the medium wide one, and vice versa 
along the x axis. 

Fig. 6 shows the total rotation measured by the 
VO system vs. the one imposed by the rotary 
stage in the rotation test. Fig. 7 illustrates the 
measured errors, which are the differences 
between the measured and the imposed total 
rotations in each position. Fig. 8 depicts the 
standard uncertainties evaluated by the Monte 
Carlo approach for the total rotation around the 
vertical axis y, while Fig. 9 shows the standard 
uncertainties evaluated for each rotational motion 



step. The uncertainties of Fig. 8 are obtained 
combining together those depicted in Fig .9. 

Measured vs. imposed total rotation 
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Fig. 6: Measured total rotation in the rotation test 



Total rotation error vs. imposed displacement 
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Fig. 7: Total rotation error in the rotation test 
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Fig. 8: Total rotation standard uncertainty around 
y axis in the rotation test 



In the rotation test, the VO system equipped with 
the telephoto lens is able to measure the rotation 
up to an angle of about 20°. This behavior is due 
to the fact that the tests are performed in a 
rectangular room: when the rotation angle is about 
0° the stereo camera optical axes are parallel to 
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the long walls of the room, while at 90° they are 
parallel to the short walls. Thus, the objects 
observed by the telephoto lenses at 0° are located 
at a distance of about 13m, but the rotation from 
0° towards 90° makes the observed objects closer 
to the cameras, up to a distance of about 2-3 m. 
This shorter distance of the observed objects 
means the number of correctly matched features 
(3D points) decreases, until the VO algorithm is 
not able to evaluate the motion. 

Rotation standard uncertainty around y axis 
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Fig. 9: Rotation standard uncertainty (of each single 
motion step) around y axis in the rotation test 

Moreover, also for small total angles the 
evaluated uncertainties (Fig. 8 and 9) are higher 
for the telephoto lens than with the other two 
lenses. These results confirm that the telephoto is 
the less suitable lens for VO in a closed 
environment as observed in [28] . 

Possible reasons of the poor performance of the 
telephoto lens are: 

• The intrinsic and extrinsic parameters 
evaluated by the Zhang approach [27] for the 
telephoto lenses are affected by wider 
uncertainties than those obtained for the other two 
focal lengths. 

• A telephoto lens has a narrower depth of 
field than a wide angle lens. 

• The number of acquired and correctly 
matched 3D points at each motion step is low due 
to the small (for a telephoto lens) distance of the 
observed objects 

Also in the rotation test, the two wide angle 
lenses behave very well. The medium wide angle 
lens has a slight advantage up to an angle of 70° . 
At high rotation values, when the distance of 
observed objects approaches its minimum values, 
both the wide angle lenses begin to exhibit higher 
uncertainties. However, the medium wide angle 
lens (with a narrower field of view) sensibly 
worsens from about 70°, while the ultra wide 
angle starts later to increase its uncertainty, from 



about 80° (see Fig. 9). This means that the 
advantage of the medium wide angle over the 
ultra wide one seems to disappear in narrow 
spaces . This effect is not present in the translation 
test, since in that case there is not a large 
reduction of the objects distance during the 
motion. 

Possible reasons of the slight discrepancy 
between the simulations performed in [4] and the 
presented experimental results could be: 

• the experimental tests are performed 
inside a closed environment. In this set-up the 
distance between the acquired objects (3D points) 
and the stereo cameras is necessarily limited by 
walls, the floor and ceil. 

• the wide difference between the 
horizontal and vertical field of view. 

5. Conclusions 

A detailed experimental comparison among 
visual odometry systems using lenses with three 
different focal lengths (an ultra wide angle, a 
medium wide angle and a telephoto lens) was 
performed and described. The experimental 
results allow to point out that in a closed 
environment the two wide angle systems 
outperform the telephoto one and that the 
difference in the field of view can lead to slight 
different performances also between the ultra 
wide angle and the medium wide angle systems . 
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Abstract 

Image quality assurance tests on medical 
ultrasound systems are often performed by 
technicians using ultrasound phantoms and 
results depend on phantom features as well as 
scanner settings and operator experience. The 
aim of the present study is the evaluation of 
variations on some features of the B-mode image 
when the ultrasound probe is handled by the 
technician during a routine quality test: 
ultrasound images of two different ultrasound 
phantom are acquired and processed from two 
transducer to evaluate measurement dispersion in 
spatial resolution, penetration depth and 
accuracy in distance measurements when probe is 
handled by the operator. Results are then 
investigated and discussed. 

1. Introduction 

Image quality assurance tests in medical 
ultrasound system are usually based on 
measurements of specific parameters such as 
spatial resolution, low contrast penetration, 
accuracy in distance measurements and local 
dynamic range [1-7]: routine tests are often 
performed by technicians using ultrasound 
phantoms 1 and results depend on phantom 
features as well as scanner settings and operator 
manual skill and experience. In scientific 
literature the operator influence on tests precision 
seems not to be systematically evaluated and in 
particular effects of probe position on 
measurements of the above parameters: in this 
work the influence of probe position is estimated 
while the transducer is handled by the technician 
during the test. 



1 An Ultrasound Phantom is a volume of material 
behaving in essentially the same manner as tissue of 
the same dimensions, with respect to absorption and 
scattering of the ultrasound radiation in question, used 
for dosimetry or for the evaluation of sonographic 
images in diagnostic sonography. 



2. Experimental set-up and methods 

The measurement setup is made of an ultrasound 
scanner Philips HD3, an array probe on the 
ultrasound phantom, and a notebook pc for data 
acquisition and processing. Tests are performed 
with two different models of array probe (convex 
and linear array) on two tissue mimicking 
phantoms 2 (see table 1 for specifications). In 
particular the study focuses on evaluation of the 
follow features: (a) accuracy in distance 
measurements (both vertical and horizontal 
directions) [3, 5], (b) high contrast spatial 
resolution (both lateral and axial resolution) [5] 
and (c) maximum depth of signal visualization [6- 
8]. 



Table 1 - Ultrasound Phantoms characteristics 



Model GAM M EX 405 GSX LE 


GAMMEX 1425 A 




(phantom A) 


(phantom B) 


Dimensions 


23.2x8.25x18.5 cm 


40.7 x 22.9 x 35.6 cm 


Weight 


Approx. 2.8 kg 


Approx. 10 kg 


Background Material 


Speed of sound 1 5 40± 1 0 m/s 


1540±10m/s 


Attenuation 


0.7±0.05 dB/cm/MHz 


0.5±0.05 dB/cm/MHz 


Pin targets 


Diameter 


0.1 mm 


0.1 mm 


Vertical 
spacing 


2 cm at 2-16 cm deep 


2 cm at 3-17 cm deep 


Horizontal 
spacing 


3 cm at 2-12 cm deep 


3 cm at 3-13 cm deep 





The ultrasound probe is applied on the phantom 
by hand and its position is manually adjusted to 
clearly display test objects in the US image, after 
that, a single image is acquired and then exported 
on a peripheral in a TIFF or DICOM format 
(lossless): for each phantom the whole procedure 
is repeated 16 times, so for a same test (see tab.2 
for settings) 16 different ultrasound image are 
provided and processed in a Matlab environment 



2 tissue-mimicking material: material in which the 
propagation velocity (speed of sound), reflecting, 
scattering and attenuating properties are similar to 
those of soft tissue for ultrasound in the frequency 
range 0,5 MHz to 15 MHz. 
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[3, 5-9] to evaluate parameters (a), (b) and (c) 
variation due to probe manual application on the 
phantom surface. 

Table 2 - Measurements settings. Test 1,2: Maximum 
Depth of Penetration, Accuracy in distance 
measurements (vertical and horizontal), High contrast 
spatial resolution (lateral and axial). Test 3, 4: 
Accuracy in distance measurements (vertical and 
horizontal), High contrast spatial resolution (lateral 
and axial). For each test and phantom 16 ultrasound 
images are acquired at a medium overall gain and a 
linear post processing. 







PROBE 1 (convex array) 






Testl 


Test 2 


Test 3 


Test 4 


Nominal 










Frequency 

range 

(MHz) 


2-5 


2-5 


2-5 


2-5 


Field of 










View 


180 


180 


110 


110 


(mm) 










Focus 










depth 


70 


70 


70 


70 


(mm) 










Dynamic 
Range 


Maximum Medium 


Maximum 


Medium 


(dB) 










Phantom 


A, B 


A, B 


A, B 


A, B 



PROBE 2 (linear array) 

Testl Test 2 Test 3 Test 4 

Nominal 



range 

(MHz) 

Field of 

View 75 75 50 50 

(mm) 

Focus 

depth 20 20 20 20 

(mm) 

Dynamic 

Range Maximum Medium Maximum Medium 

(dB) 

Phantom AJ$ AJ$ AJ$ A, B 



To better understand settings in table 2 and before 
showing measurements and related results, some 
definitions are reported below [5-8]: 

• Nominal Frequency (of a transducer): acoustic 
working frequency of a transducer as quoted 
by the designer or manufacturer. In our study 
the ultrasound scanner does not provide a 
single value for nominal frequency but a range, 
i.e. 2-5 MHz for the convex array probe, 6- 
9MHz for the linear array probe. 

• Field of View (FoV): area in the ultrasonic 
scan plane from which ultrasound information 
is acquired to produce one image frame. Since 
for each probe the size in pixel of the 
diagnostic image is quite constant a FoV 
variation produces a different scale factor (mm 



to pixel ratio), influencing other image 
characteristics, i.e. spatial resolution. 

• Focus Depth: nominal depth of the 
transmission focus in the ultrasound image. 
Transmission foci are usually marked on one 
side of the diagnostic image. 

• Dynamic range (global): ratio of the maximum 
to the minimum echo -signal amplitude, even 
with changes of settings, that a scanner can 
process without distortion of the output signal. 

• Local Dynamic Range: ratio, expressed in 
decibels, of the minimum echo amplitude that 
yields the maximum grey level in the digitized 
image to that of the echo that yields the lowest 
grey level at the same location in the image 
and the same settings. Local Dynamic Range 
influences image contrast and is setting 
dependent (i.e. through the dynamic range 
control on the ultrasound scanner). 

• Gray Scale Mapping Function (GSMF): 
relationship between echo amplitudes and gray 
levels on the image. From GSMF the Local 
Dynamic Range can be estimated. 

Settings have been chosen to display the 
maximum number of test object at the central 
frequency of the ultrasound probe (working 
frequency). Moreover all images are acquired 
within ±2° tilt angle, at medium overall gain and 
linear post processing while persistence, edge 
enhancement and other image processing are set 
to minimum. 



2.1 Accuracy in distance measurements 

Wires embedded within the ultrasound phantoms 
are imaged with the sensitivity adjusted to make 
the displayed echoes as sharp as possible: 
horizontal and vertical distances are measured 
"from peak to peak" of wires sections and 
compared with their nominal values by software 
[3]. In particular, after the ultrasound image is 
acquired and processed in a workstation for 
setting the scale factor in the diagnostic area, the 
operator is asked to choice test object pairs: it is 
not necessary centering each target because the 
software developed automatically calculates its 
barycentric coordinates within a Region of 
Interest (ROI) and uses it for measuring the 
distance among each pair of wires, so errors due 
to operator's visual acuity can be neglected. For 
each pair the measurement distance relative error 
e%=|d r k-dk|/d r k x 100 between the nominal distance 
d rk and the measured value d k in the image is 
calculated for both horizontal and vertical 
distances. 
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2.2 High Contrast Spatial resolution 

High-contrast resolution characteristics can be 
obtained by measuring dimensions of the point- 
spread function (PSF), which is the characteristic 
response of the imaging system to a high-contrast 
point target: PSF can be produced by imaging 
high reflection targets that are smaller than one 
wavelength, as nylon wires section in PSF 
phantoms of table 1 . For ultrasound scanners the 
PSF is not singular, nor isotropic and has different 
axial and lateral size in the Field of View (FoV), 
depending on distance from the transducer 
emitting surface, so high contrast spatial 
resolution also changes with position and depth in 
the image. Thus, many different measurements of 
the axial and lateral diameter of displayed wires 
sections at different positions and depths have 
been performed by software for each settings in 
table 2 to obtain representative values of the 
system's axial and lateral resolution. In particular 
for each wire section k (point target) in the US 
image a threshold algorithm at Full Width Half 
Maximum is applied to evaluate both its 
maximum width PSFx max (k) and length 
PSFz max (k): PSFx max (k) and PSFz max (k) values are 
assigned as high contrast axial and lateral 
resolution respectively 



2.3 Maximum depth of signal visualization 

Maximum Depth of Signal Visualization DSV max 
is the maximum depth at which echo signals from 
scattering within a tissue-mimicking phantom can 
be detected [6-8]. In particular, as shown in fig. 1, 
for each image DSV max is calculated from a depth 
profile (mean gray level vs depth) of gray levels 
within a ROI DP of 30 pixel width, by means of a 
threshold at 2dB [8] above the mean value 
displayed in a ROI n of 10x10 pixels placed at the 
end of the Field of View (which is associated to 
the electronic noise level), applied to the US 
image. The +2dB threshold needs the relationship 
between gray level in the image and acoustical dB 
[8], so a Gray Scale Mapping Function is 
estimated by means of the procedure in [9] for 
each probe and dynamic range setting in table 2. 
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Fig. 1: Maximum Depth of Signal Visualization: 
DSV m ax is evaluated from the depth profile of gray 
levels in ROI DP at +2dB over the noise level. In 
ROI n the mean gray level estimates the noise level. 



3. Results and Discussion 

For each of the above parameters results are 
classified for probe and phantom model. Since all 
measurements are done by means of an image 
analysis software that minimize errors due to 
operator's visual acuity and subjective judgment, 
principal causes of results dispersion can be 
related to small image variations due to manual 
probe position on the phantom surface. Therefore 
in our study influences of ultrasound transducer 
position on quality assurance measurements are 
estimated by means of expanded uncertainties, 
evaluated as measurement repeatability [10] with 
a coverage factor t= 2.13 at 95 percent confidence 
level. 



3.1 Accuracy in distance measurements 

Accuracy in distance measurements has been 
evaluated for all the settings in table 2. In figure 2 
results for the convex prove are shown: for both 
vertical and horizontal measurements no 
significant differences between tests have been 
observed, while for horizontal distances (fig. 2b 
and fig. 2d) discrepancies between phantoms can 
be noticed. Moreover, to evaluate test influences 
on measurement repeatability a Fisher test (F-test) 
has been performed for results at vertical- 8 0mm 
and horizontal- 6 0mm: while in vertical distance 
measurements no significant differences can be 
confirmed among variances, in horizontal ones 
homoscedasticity should be rejected (for both of 
phantoms), so we can suppose a variation in 
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measurement repeatability related to test settings. 
Measurements uncertainty ranges from: 0.3 to 0.7 
percent for vertical distances, 0.2 to 1.3 percent 
for vertical distances. 




£ 2 



20 40 60 80 100 120 140 
nominal distance - vertical (mm) 
(a) 




30 60 90 

nominal distance - horizontal (mm) 
(b) 




20 40 60 80 100 120 
nominal distance - vertical (mm) 
(c) 



140 




30 60 90 

nominal distance - horizontal (mm) 
(d) 



'Test 1 ' 



'Test 2 



Test 3 ' 



»Test 4 



Fig. 2: Distance measurements for convex array 
probe. Phantom A: (a) vertical distances, (b) 
horizontal distances. Phantom B: (c) vertical 
distances, (d) horizontal distances. 

In figure 3 results are shown for the linear array 
probe and no significant differences are observed 



for most of the measurements, nevertheless the F- 
test indicates no equality of variances on vertical 
(20 mm, phantom A) and horizontal distances 
(30mm, both of phantoms). Measurements 
uncertainty ranges from: 0.1 to 5.5 percent for 
vertical distances, 0.1 to 0.6 percent for 
horizontal distances. 
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Fig. 3: Distance measurements for linear array 
probe. Phantom A: (a) vertical distances, (b) 
horizontal distances. Phantom B: (c) vertical 
distances, (d) horizontal distances. 

For convex array transducer differences between 
test in horizontal distances (fig.2b and fig.2d) can 
be due to phantom tissue mimicking material (e.g. 
phantom A attenuation is 1.4 times greater than in 
phantom B so pin target imaging could be more 
smooth and less visible over the background 
speckle) and mostly to Field of View variation, 
because it influences mm to pixel ratio (scale 
factor) and the distortion of imaged nylon wires 
sections. In fact at lower FoVs, wires at 60mm are 
selected near image boundaries (at 20-3 0mm 
depth), where distortion artifacts are more intense 
and become more noticeable for higher contrasts, 
i.e. depending on the Dynamic Range settings. All 
of the causes above likely influence the 
measurement repeatability, enhancing slight 
differences in acquired images and their 



59 



IX Congresso MMT 



processing by the analysis software. Similar 
considerations can be supposed for the linear 
array transducer, where the presence of high noise 
levels adds dispersion to results. 



3.2 High Contrast Spatial resolution 

Axial and lateral resolution at FWHM have been 
evaluated for all the settings in table 2. In figure 4 
no significant variations between tests have been 
observed for most of axial and lateral 
measurements of the convex array probe, 
nevertheless some discrepancies between test can 
be noticed in phantom A. A Fisher test has been 
performed for results at 100 mm (both for axial 
and lateral resolution) and indicates that for 
results in fig.4c and fig.4d (phantom B) 
homoscedasticity should be rejected. 
Measurements uncertainty ranges from: 0.1 mm 
to 1.0 mm for axial resolution, 2.9 mm to 7.9 mm 
for lateral resolution. 

In figure 5 spatial measurements are shown for 
the linear array probe and no significant 
differences are observed for most of the 
measurements, nevertheless the F-test indicates 
no equality of variances on both axial and lateral 
resolution values obtained with the phantom B. 
Measurements uncertainty ranges from: 0.1 mm 
to 0.4 mm for axial resolution, 0.1 mm to 1.4 mm 
for lateral resolution. 

Differences between tests can be due to noise 
level, phantom tissue mimicking material and 
Field of View variation, because of it influence 
on mm to pixel ratio, and so both on scale factor 
and on pixelation error. As in distance 
evaluations, measurement repeatability likely 
depends on all of the causes above, because of the 
consequent enhancement of differences among 
acquired images. 
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Fig. 4: High contrast spatial resolution 
measurements for convex array probe. Phantom A: 
(a) axial resolution (b) lateral resolution. Phantom 
B: (a) axial resolution (b) lateral resolution 
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DSV max is more settings dependent while 
discrepancies between phantom are not 
significant: we suggest it can be due to excessive 
noise increasing with contrast (Dynamic Range 
variation). Variations in measurements 
repeatability could be due to phantom attenuation, 
e.g. higher attenuation provides less dependence 
on image tilting and consequent artifacts. 
180 




phantom A phantom B 

■ Test 1 ■ Test 2 

Fig. 6: Maximum depth of signal visualization for 
convex array probe 



♦ Testl ■ Test 2 Test 3 X Test 4 

Fig. 5: High contrast spatial resolution 
measurements for linear array probe. Phantom A: 
(a) axial resolution (b) lateral resolution. Phantom 
B: (a) axial resolution (b) lateral resolution 

3.3 Maximum depth of signal visualization 

Maximum depth of signal visualization DSV max 
has been evaluated for settings test 1 and test 2 in 
table 2. In figure 6 a significant discrepancy 
between phantom is observed for both of tests: 
higher depth values are related to phantom B. The 
F-test indicates that homoscedasticity should be 
rejected only for phantom B results. 
Measurements uncertainty ranges from 2 mm to 
5 mm. 

In figure 7 a significant difference in results 
between tests and phantom can be noticed: higher 
depth values are consistent to each other and 
related to settings of test 2 for both of the 
phantoms. Results from F-test denotes that 
variance should be considered non-homogeneous 
between tests for both of phantoms. 
Measurements uncertainty ranges from 1 mm to 
10 mm. 

In convex array analysis DSV max seems to depend 
mostly on phantom model than test settings, it is 
likely due to different attenuation coefficient 
(phantom B attenuates less than phantom A). On 
the other hand, in linear array examinations 
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phantom A phantom B 
■ Test 1 ■ Test 2 

Fig. 7: Maximum depth of signal visualization for 
linear array probe 

4. Conclusion 

In the present study variations on some features of 
the B-mode image have been investigated when 
the ultrasound probe is handled by the technician 
during a routine quality test: ultrasound images of 
two different ultrasound phantom are acquired 
from two ultrasound transducer and processed by 
software to evaluate measurement dispersion in 
high contrast spatial resolution, penetration depth 
and accuracy in distance measurements. 
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Measurement uncertainties depend strongly on 
probe model and parameter investigated: they 
ranges from 0.1 to 5.5 percent in distance 
measurements error, from 0.1 to 7.9 mm in high 
contrast resolution and from 2 to 10 mm in 
maximum depth of signal visualization. Although 
numerical results are limited to the two examined 
probes some general consideration can be done: 
(a) measurements depend on settings as well on 
phantoms features, probes and parameters 
investigated (b) Field of View settings seems to 
be more important for measurement repeatability 
than Dynamic Range, (c) accuracy in distance 
measurements seems to be less dependent on 
settings and transducer position. However use of 
support stands and clamps is highly recommended 
during routine quality assurance tests in order to 
reduce uncertainty in ultrasound system 
performances estimation. 

Other measurements are going to be collected on 
different ultrasound systems to perform an in- 
depth analysis of uncertainties in diagnostic 
ultrasound quality assurance by means of 
ultrasound phantoms and image analysis software. 
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Abstract 

This work describes the design of two contactless 
temperature measurement methods based on total 
radiance and two-colour pyrometry. The methods 
have been conceived to measure the temperature of 
a small brass coated steel wire during wire 
drawing. Usual contact sensors cannot be applied 
due to the wire movement and excessively large 
insertion errors in this critical condition. The 
pyrometers optical layouts have been analysed by 
means of numerical simulations in order to 
evidence their sensitivity to the wire oscillations. 
Performances of different infrared sensors have 
been compared on the basis of the achieved 
measurement uncertainty simulating background 
temperature variation, slope of the wire's 
emissivity and the effect of the atmosphere 
absorption. 

1. Introduction 

Infrared thermometers allow temperature 
measurement without any contact between the 
sensor and an observed object. This is a clear 
advantage in case the object size is so small that 
even with the smallest contact thermometers (e.g. 
miniature thermocouples or RTDs) the loading 
effect due to the sensor cannot be accepted. 
Anyway, contactless temperature measurement is 
not an easy task since achievable accuracy relies 
directly on the knowledge of many quantities such 
as the object emissivity, the background 
temperature, the medium spectral transmissivity 
and on the condition that the emitting source fills 
completely the sensor Field of View (FOV) [1-5]. 
Emissivity is roughly known because depends on 
temperature and might change because of the local 
state of oxidation, cleanliness or roughness [2] of 
the observed surface. These criticalities affects the 
analysed case study, i.e. wire drawing, since the 
manufacturing process requires soaps and 
lubricant oils and the wire emissivity varies 
because of the material inhomogeneity or local 
oxidation, therefore the usage of common total 
radiance pyrometers can lead to large measurement 
errors. Moreover, the wire size (between 0.1 to 1 



mm) is generally smaller than the measurement 
spot of any standard pyrometer or infrared camera 
[6]. This leads to measurement errors because 
beside the wire' emission some radiance from the 
background scene is measured as well, and 
averaged with the source one. Finally, significant 
transversal oscillations of the wire are present 
during the wire drawing and brass coating 
deposition process. In case of wire oscillations, 
radiative coupling between the source and the 
sensor changes with the wire position, causing the 
fluctuation of the measured temperature or in 
general, a SNR (Signal to Noise Ratio) worsening. 

Thus, two contactless temperature measurement 
methods have been developed to overcome the 
drawbacks of the standard pyrometers and 
minimize the achievable uncertainty in measuring 
the temperature of small wires with large 
oscillations and emissivity changes. The expected 
improvements with respect to off the shelf 
instruments derive from purposely designed 
instruments optical layouts and the optimal choice 
of the detectors. 

The paper is organized as follows: Section 2 
summarizes the modelling of the two pyrometers 
whereas Section 3 describes the conceived optical 
layouts. Section 4 deals with the uncertainty and 
sensitivity analyses performed for the proposed 
instruments whereas Section 5 summarizes the 
numerical results. The obtained results are 
discussed in Section 6 and the manuscript is 
eventually concluded in Section 7. 

2. Backgrounds and methods 

According to the Planck law, the heat flux for unit 
wavelength at a given temperature is defined as: 

E nX = % Eq.l 

where k is the wavelength whereas C/ and C2 are 
known constants. In order to account for the 
emissivity of the source and the efficiency of the 
radiative exchange between two surfaces i and j, 
the radiative conductor (GR) is defined as follows: 
GR t =e l A l GB t Eg. 2 



63 



IX Congresso MMT 



where A% is the emitting surface area, Si is the 
surface emissivity and GBy is the Gebhart factor 
that expresses the fraction of power radiated from 
surface i and absorbed by surface j. The latter 
factor accounts for the multiple reflections of the 
wire radiance reaching the detector. 

Finally, the detector output O for a source at a 
given temperature Tis: 

O(T) = £ maX R(X ) GR E nk (T) dX Eq. 3 

where k max and k m i n are the extremes of the 
detector wavelengths range, R(X) the detector 
sensitivity, E n x (T) the source spectral radiance. 

Conversely, the two-colour pyrometer 
measurement method is based on the ratio between 
the heat fluxes at two different wavelengths; 
considering the two sensors outputs one can write: 



o, 



GRj R 1 (A 1 )s(A 1 )A 2 



GR 2 R 2 (A 2 )e(A 2 )A 1 5 



c 2 



-1 



Eq. • 



The measured temperature can be derived from 
Equation 4 under the assumption that the source 
emissivity is a smooth function, and so, similar at 
the two adjacent wavelengths: 

f I 1^ 



^A 2 A 1 



C 7 



r E nM 



In 



F A 



Eq. 5 



Moreover, Equation 5 is valid under the 
hypothesis that the exponential term in Equation 4 
is at least one order magnitude larger than one. The 
latter assumption is fully satisfied within the 
investigated wavelength range, i.e. between 0.5 
and 5 jam, and at the analysed temperatures (from 
200°C to 800°C). 

3. Optical layouts design 

The radiative models of the pyrometers have been 
developed with ESATAN-TMS software. The 
model includes the infrared sensor active area, the 
mirrors, the wire and an enveloping box. Figure 1 
shows the conceptual optical layouts. The first 
configuration for the total radiance pyrometer 
exploits one planar diffusive mirror to reflect the 
radiance emitted by the wire. The mirror is 
replaced in the two-colour pyrometer by parabolic 
reflectors. The internal surfaces of the optical 
cavity (mirrors and envelope box) aim to reflect the 
wire radiance that does not directly reach the 
sensor and therefore, materials with high infrared 
reflectivitv have been used. 



Planar Parabolic 
^^^^H Mirror ^^^k Mirror 

Sensor > ,^ 



Fig. 1: Pyrometers optical layouts; (left) total 
radiance configuration, (right) optical configuration 
for the two-colour pyrometer. 

The optical configurations have to be compared 
on the basis of the radiative efficiency and 
sensitivity to the wire oscillations. 

Thus, GRs between wire and detector have been 
computed in nine different positions, as evidenced 
in Figure 2. 



Wire' positions 





Detector 

»8 ■ 







Fig. 2: Source position with respect to the detector. 

The maximum distance from the reference 
position (position "0" in Figure 2) is 0.5 mm, value 
obtained from the measurement of the wire 
oscillation amplitudes in preliminary tests. The 
optical properties (emissivity, transmissivity and 
specular or diffusive reflectivities) of the materials 
are reported in Table 1 . In particular, the coating 
for the total radiance pyrometer is sand blasted 
aluminium whereas gold is used in the two-colour 
configuration. The meshes of the optical elements 
in each configuration have been defined by 
preliminary analyses aiming to identify the 
maximum size for each component granting 
stability in the ray tracing analysis. 



Material 


£ 


T 


p 

specular 


P 

diffusive 


Brass 


0.4 


0 


0 


0.6 


Sand blasted 
aluminium 


0.15 


0 


0 


0.85 


Gold IR 


0.03 


0 


0.9 


0.07 



Table 1 Materials optical properties 

Table 2 summarizes geometry and mesh for the 
analysed optical layouts. 
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Component 


Mesh 
size 


Model 

Size 


Wire 


30x30 


32 mm, diameter 0.8 
mm 


Envelope Box 


22x22x22 


30x30x26 mm 3 


Planar Mirror 


15x15 


30x30 mm 2 


Parabolic 


40x40 


Diameter 30 mm, 


Mirror 1 


Focal length 36 mm; 


Parabolic 


40x40 


Diameter 30 mm, 


Mirror 2 


Focal length 26 mm; 


Detectors 
Active area 


lxl 


lxl mm 2 



Table 2 Models geometry and mesh size. 

The radiative analyses have been performed in 
each position of Figure 2 with the optical 
properties of Table 1. Figures 3 shows the 
computed GRs. Numerical simulations evidenced 
that for the total radiance pyrometer the planar 
mirror configuration with sand blasted aluminium 
coating allows reducing the sensitivity to the wire 
oscillations. In fact, the average GR computed 
among the different wire positions is about 
2.3 1 0" 8 m 2 , with a standard deviation 0.85% of the 
average value. 

2.35E-08 





Wire position 

Fig. 3: Computed GRs between wire and detector in 
the total radiance pyrometer (top) and in the two- 
colour configuration (bottom). 

Conversely, in the two-colour pyrometer, the 
configuration with two parabolic mirrors and gold 
coating achieves the highest GR, providing a 
maximum of about 7.54 10" 8 m 2 but at the price of 



a larger variability (about 4% with respect to the 
average value). 

4. Temperature uncertainty 
evaluation 

In order to evaluate the expected uncertainty for 
the measured temperature, the first contribution 
comes from the detector noise [7], i.e. the photon 
noise. Typically, datasheets of commercial 
detectors provide the responses and the 
sensitivities of the detectors in terms of spectral 
detectivity D*(o): 



D(a) 



Eg. (6) 



NEP(a) 

where Ad is the detector area, Af is the sensor 
frequency bandwidth and NEP is the detector noise 
equivalent power. Both the D*(a) and the NEP 
depend on the wavenumber o. In our case, once Af 
has been set (i.e. 1 Hz), the NEP can be found 
reversing Equation 6. The ratio between the 
measured power and the NEP gives the achievable 
SNR. 

Beside the uncertainty introduced by the photon 
noise, the radiance of the background leads to an 
overestimation of the source radiance. This effect 
has been evaluated in both pyrometers. Moreover, 
the variability of the total or spectral wire 
emissivity is expected to be of paramount 
importance to define the final measurement 
accuracy. Thus, some sensitivity analyses have 
been performed to identify the variability of the 
measured temperature due to a global emissivity 
change in the total radiance pyrometer, and a 
spectral emissivity slope for the two-colour 
configuration. 

Three photoconductive infrared detectors have 
been analysed, i.e. Silicon (Si), Lead Sulfide (PbS) 
and Lead Selenide (PbSe), covering different 
wavelength ranges, responsivities and 
detectivities. The detectors characteristics are 
summarized in Table 3. 



Type 




Si 


PbS 


PbSe 


Case 




TO- 18 


TO- 18 


TO- 18 


A d 


mm 2 


1 


1 


1 


Range 


[xm 


0.4-1.2 


1-2.8 


K4.5 


D* 


cm Hz 0 5 W 1 


4.4 10 12 


1 10 11 


5 10 10 


NEP 


W 


2.3 10" 14 


1 10" 12 


2 lO" 12 



Table 3 Characteristics of the analysed detectors. 

The uncertainty and sensitivity analyses have 
been performed for the total radiance pyrometer 
with the following parameters: 
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• The nominal wire temperature has been set 
to 200, 400, 600 and 800 °C; 

• Atmosphere spectral absorption at the 
different temperatures of the wire has been 
modelled [8]; 

• Worst case background emission has been 
simulated, with unitary emissivity and 
temperature (T BK ) set at 30, 35 and 40 °C; 

• A normal distribution of the total 
emissivity of the wire has been considered 
to account for the emissivity variation 
because of oils, soaps or oxidation in the 
manufacturing process. The mean 
emissivity has been set to 0.4 and the 
standard deviation to 10% of the average 
value. 

The analyses in the two-colour pyrometer have 
been performed with the following additional 
parameters: 

• The transmittance curves of commercially 
available infrared filters have been added 
to compute the measured power; 

• A variable wire emissivity difference 
between the two reference wavelengths, 
up to 5%, has been simulated to take into 
account a possible variation of the source 
spectral emission; the simulated value is in 
agreement with literature data for brass, 
i.e. 20% change over 2.7 jam, as shown in 
reference [9]. 

5. Analysis results 

Hereafter the results of the performed simulations 
are described for both the analysed pyrometers. 

5.1. Total radiance pyrometer 

The SNRs for the selected sensors at the different 
temperatures are reported in Table 4 whereas Table 
5 provides the expected temperature variability 
because of the wire oscillations. 



T[K] Si PbS PbSe 



473 


3.8 


2.6 10 5 


3.1 10 6 


673 


2.1 10 4 


1.1 10 7 


3 10 7 


873 


2.7 10 6 


9.1 10 7 


1.3 10 8 


1073 


6.1 10 7 


3.8 10 8 


3.4 10 8 



Table 4 SNRs at different temperatures of the 
analysed detectors. 

The temperature error introduced by the 
background radiance is summarized in Figure 5 
whereas, the temperature uncertainty related to the 
wire emissivity variability is highlighted in Table 
6. In the latter case, the temperature uncertainty has 



been derived using the numerical propagation as 
evidenced in reference [10]. 



T 


<TSi 


<TPbS 


OPbSe 


473 


n.d. 


0.047 


0.12 


673 


0.067 


0.13 


0.13 


873 


0.087 


0.16 


0.16 


1073 


0.11 


0.33 


0.22 



Table 5 Variability of the measured temperature 
due to wire oscillation. Temperatures are in kelvin. 
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303 308 313 
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Fig. 5 Expected temperature error due to 
background emission with the Si (top), PbS 
(middle) and PbSe (bottom) detectors. 



T[K] 


Si 


PbS 


PbSe 


473 




0.69 


1.06 


673 


0.49 


1.14 


1.69 


873 


0.61 


1.38 


1.96 


1073 


0.72 


1.47 


2.10 



Table 6 Relative uncertainty of the measured 
temperature because of the wire' emissivity 
distribution. The uncertainty is expressed as 
percentage. 
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l.E-04 




I.E-08 1 1 ' 

473 673 873 1073 

T[K] 

Fig. 6 Infrared filters comparison for Si (top), PbS 
(middle) and PbSe (bottom) detectors. 



Filter 


CWLs 






T[K] 




pair 


[nm] 


473 


673 


873 


1073 


sO-sl 


730-766 


0.36 


0.49 


0.62 


0.74 


sl-s2 


766-800 


0.42 


0.59 


0.75 


0.91 


s2-s3 


800-830 


0.59 


0.88 


1.19 


1.53 


s3-s4 


830-852 


0.74 


1.03 


1.30 


1.56 


s4-s5 


852-855 


70.60 


24.38 


14.43 


11.53 


s5-s6 


855-880 


0.70 


0.98 


1.25 


1.51 


s6-s7 


880-905 


0.95 


1.50 


2.17 


3.02 


s7-s8 


905-940 


0.62 


0.90 


1.19 


1.49 


s8-s9 


940-1064 


0.26 


0.42 


0.63 


0.93 



Table 6 Relative temperature error due to wire 
oscillation for different filter pairs, Si detector. 
Error is expressed as percentage. 



Filter 


CWLs 




T[K] 




pair 


[jim] 


473 


673 


873 


1073 


tO-tl 


1.25-1.3 


0.75 


1.04 


1.32 


1.58 


tl-t2 


1.3-1.35 


0.81 


1.12 


1.42 


1.70 


t2-t3 


1.35-1.4 


0.86 


1.20 


1.51 


1.81 


t3-t4 


1.4-1.45 


0.92 


1.28 


1.61 


1.93 


t4-t5 


1.45-1.5 


0.99 


1.36 


1.72 


2.05 


t5-t6 


1.5-1.55 


1.05 


1.45 


1.82 


2.18 


t6-t7 


1.55-1.6 


1.12 


1.54 


1.93 


2.30 


t7-t8 


1.6-1.65 


1.18 


1.63 


2.04 


2.43 



Table 7 Relative temperature error due to wire 
oscillation for different filter pairs, PbS detector. 
Error is expressed as percentage. 



Filters 


CWLs 






T[K] 




pair 


[jim] 


473 


673 


873 


1073 


uO-ul 


2.7-2.95 


0.77 


1.11 


1.44 


1.77 


ul-u2 


2.95-3.46 


0.38 


0.48 


0.56 


0.63 


u2-u3 


3.46-3.6 


1.98 


2.95 


3.79 


4.56 


u3-u4 


3.6-4.26 


0.57 


0.76 


0.95 


1.13 


u4-u5 


4.26-4.67 


1.37 


2.44 


3.69 


5.24 



Table 8 Relative temperature error due to wire 
oscillation for different filter pairs, PbSe detector. 
Error is expressed as percentage. 



5.1. Two-colour pyrometer 

Figure 6 compares the computed relative 
uncertainty in temperature measurement for some 
commercially available IR filters. The uncertainty 
has been computed propagating the uncertainty 
contributions of Equation 5. The uncertainty about 
the heat flux measured by each detector has been 
evaluated as the ratio between the detector' NEP 
and the measured radiance. The central 
wavelengths (CWLs) of the analysed filters range 
between 730 and 1064 nm for the Si detector, from 
1250 to 1650 nm for the PbS and within 
2700-4670 nm for the PbSe. 

The same filters have been used to highlight the 
pyrometer sensitivity against the wire oscillations. 
Results about the expected temperature variability 
are summarized in Table 7, Table 8 and Table 9. 

The relative error caused by variations of the wire 
spectral emissivity within the IR band filters is 
evidenced in Figure 7. The error has been 
computed only considering the infrared filters 
selected for each detector (see the filters pairs in 
green bold in Tables 6, 7 and 8). Finally, the error 
due to background radiation is summarized in 
Figure 8. 
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Emissivity slope 

Fig. 7 Temperature error due to spectral 
emissivity slope between the selected filters pairs 
for the Si (top), PbS (middle) and PbSe (bottom) 
detectors. Red dashed lines evidence the expected 
emissivity differences for the selected filters pairs 
for each detector. 




Fig. 8 Temperature relative error due to the 
background emission for the Si (top), PbS (middle) 
and PbSe (bottom) detectors. 

6. Discussion 

The computed SNRs for the analysed detectors in 
the total radiance pyrometer evidence that the best 
results are achieved with the PbSe but in general, 
high SNRs can be obtained as well for the other 
detectors. The Si provides a high SNR but only for 



temperatures larger than 473 K (see Table 4). The 
result is somehow expected since increasing the 
temperature, the Plank curves shift towards the 
lower wavelengths, i.e. below the detector cut-off 
wavelength. 

The wire oscillation in the total radiance 
configuration leads to small temperature errors for 
all detectors, as shown in Table 5. The maximum 
error is less than 0.5 K at the highest nominal wire 
temperature, i.e. 1073 K. This result is expected as 
well, since the optical layout was conceived to be 
inherently insensitive to wire oscillations, as can be 
verified by looking to the GR plots in Figure 3. 

The effect of the background emission on the 
measured temperatures is shown in Figure 5: the 
largest temperature error is obtained at the lowest 
wire temperatures for all the detectors. Moreover, 
the computed temperature error is negligible for 
the Si and PbS detectors, which show maximum 
relative errors of about 1.2 10" 8 and 0.25%. The 
error rises up to about 2.4% in case of the PbSe 
detector. This result is reasonable since the PbSe 
has a spectral range shifted to the high wavelengths 
(as shown in Fig. 4), therefore it is more sensitive 
to the background disturbance that deriving from 
low temperatures (in comparison to the wire) has 
most of its radiance in the in the mid infrared. 

The simulation about the statistical variation of 
the wire emissivity highlighted a significant 
dependence of the expected temperature 
uncertainty on the total emissivity changes (as 
shown in Table 6); in fact, the temperature 
uncertainty varies between 7.7 K, maximum value 
for the Si detector, and 22.5 K, obtained with the 
PbSe. Moreover, for all the considered detectors, 
the uncertainty increases with the increase of the 
wire temperatures. 

The analysis of some commercially available IR 
filters for the two-colour pyrometer leaded to the 
selection of the filter pairs s8-s9 for the Si detector, 
t7-t8 for the PbS and ul-u2 for the PbSe. These 
allow achieving the best SNRs (as shown in Figure 
6) and minimizing the sensitivity to the wire 
oscillations (as evidenced in Table 6, Table 7 and 
Table 8). Starting from the selected filters, the 
sensitivities about the spectral emissivity variation 
and background emission have been investigated. 
The positive and negative slopes of the wire 
spectral emissivity leads to a temperature error that 
increases with the wire temperature increase. The 
maximum computed relative error has been found 
to be less than 0.26%, 5% and 4.2% for the Si, PbS 
and PbSe, respectively. These are highlighted by 
the red dashed lines shown in Figure 7. 

The background emission is negligible for the Si 
and PbS but provides a maximum relative error of 
about 1.5%, at 200°C, for the PbSe. This error, 
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however, could be reduced if the temperature of the 
envelope box was stabilized, e.g. by using a 
thermoelectric cooler. 

7. Conclusions 

Two methods for the contactless temperature 
measurement of a thin wire during brass coating 
and drawing have been developed and optimized to 
overcome the shortcomings inherent to this 
measurement case i.e. the small wire size, 
emissivity changes, its velocity along the wire axis 
and transversal oscillations. The two methods, 
based on a total radiance and two-colour 
pyrometers, exploit commercially available 
infrared detectors and optical components. The 
instruments sensitivity to the wire oscillations, 
background emission and emissivity variation has 
been analysed through numerical modelling. 

In the total radiance pyrometer, the PbSe detector 
evidences the best performances within the 
investigated temperature range, i.e. from 200 °C to 
800 °C, providing the largest SNR and a 
temperature uncertainty between 1.8% and 2.5% 
considering the effect of emissivity variance. 

In the two-colour pyrometer, the best results are 
obtained with the PbSe detector as well. The latter 
provides the highest SNR and the lowest 
measurement uncertainty, about 0.63%, using 
infrared filters with band centres positioned 
between about 3 and 3.5 jam. In both cases, the 
optics enclosure has been assumed to be thermally 
stabilized to avoid the large errors that have been 
evidenced for the background emission. 

References 

[1] Quinghai Ren, "Design of a dual channel 
thermometer for true temperature 
measurement", International Journal of 
Infrared and Millimeters Waves, 1996, pp. 
1969-1975, Vol. 17, No. 11,doi: 
10.1007/BF02069470. 

[2] Michael F. Modest, "Radiative Heat Transfer 
Third Edition", Academic Press, 2013, ISBN 
978-0-12-386944-9. 

[3] Tairan Fu, Xiaofang Cheng, Xueliang Fan and 
Jinlei Ding, "The analysis of optimization 
criteria for multi-band pyrometry", Metrologia, 
2004, pp. 305-313, Vol. 41 , doi: 10.1 088/0026- 
1394/41/4/012. 

[4] Krzysztof Chrzanowski and Marek Szulim, 
"Measure of the influence of detector noise on 
temperature-measurement accuracy for 
multiband infrared systems", Applied Optics, 
1998, pp. 5051-5057, Vol. 37, No. 22, doi: 
10.1364/AO.37.005051. 

[5] Krzysztof Chrzanowski and Marek Szulim, 
"Errors of temperature measurement with 
multiband infrared systems", Applied Optics, 



1999, pp. 1998-2006, Vol. 39, No. 10, doi: 
10.1364/AO.38.001998. 

[6] Diego Scaccabarozzi and Bortolino Saggin, 
"About the dynamic characterization of micro- 
bolometric infrared cameras", Sensors and 
Actuators A : Physical, July 2014, doi: 
10.1016/j.sna.2014.06.018. 

[7] Davis, Sumner P., Mark C. Abrams, and 
James W. Brault, "Fourier transform 
spectrometry", Academic Press, 2001, ISBN: 
978-0-12-042510-5. 

[8] Rothman, Laurence S., et al. "The HITRAN 
molecular database: Editions of 1991 and 
1992." Journal of Quantitative Spectroscopy 
and Radiative Transfer 48.5 (1992): 469-507. 

[9] S0nnik Clausen, Axel Morgenstjerne and Ole 
Rathmann, "Measurement of surface 
temperature and emissivity by a 
multitemperature method for Fourier- 
transform infrared spectrometers", Applied 
Optics, 1999, pp. 5683-5691, 1996, Vol. 35, 
No. 28, doi: 10.1364/AO.35.005683. 

[10] JCGM 100:2008 Evaluation of measurement 
data - Guide to the expression of uncertainty 
in measurement (GUM). 



69 



IX Congresso MMT 



Measurement of the heat removed by devices for 
skin tags treatment 

Marco Tarabini*, Bortolino Saggin* 

* Politecnico di Milano, Dipartimento di Meccanica, Via Previati 1/C, 23900 Lecco 



Abstract 

This paper describes a method for the 
comparison of the cooling capability of over-the- 
counter devices for skin tags removal. The method 
is based on the principle of the heat flux-meter: 
the accuracy of temperature and heat flux 
measurements has been experimentally evaluated 
(lower than 6% for heat flux measurements and 
lower than 0.3 °C for temperature 
measurements). The instrument allowed 
comparing the thermal efficiency of nineteen 
Dimethyl ether and propane dispensers using 
different nozzles. 

1. Introduction 

Dermatological cryotherapy is used to treat 
different skin problems, such as tags, verrucas, 
molluscum contagiosum, seborrheic and actinic 
keratosis. Cryotherapy, typically performed with 
liquid nitrogen, consists of freezing a small skin 
volume for a short time period. Up to ten years 
ago, it was not clear how the low temperature 
solved the skin problem [1]. The mechanisms of 
verrucas removal are nowadays clear and have 
been clearly explained by Nguyen and Burkhart 
[2]; when the tissue temperatures falls below -5°C 
to -15°C, extracellular ice crystals begin to form. 
The crystals mechanically disrupts cellular 
membranes and disturb fluid homeostasis, 
resulting in cellular dehydration and cell death. 
Vascular changes associated with circulatory 
failure become evident between the formation of 
extracellular and intracellular ice crystals, as 
tissue temperatures fall below -15°C. 
Microthrombi form within damaged vessels 
leading to ischemic necrosis. 

The most effective cooling method is the topical 
application on lesions that one wishes to destroy 
of a cotton wool swab saturated by immersion in 
liquid nitrogen or other gases with low boiling 
temperatures. Nowadays, several over-the-counter 
(OTC) warts-freezing therapies are available. 
These devices use Dimethyl Ether and Propane 
instead of the liquid nitrogen and their efficiency 
is currently assessed by comparing their effects on 
patients, as in [3]. A very simple comparison 



between over-the-counter wart preparations and 
liquid nitrogen is presented in reference [4]. The 
comparison was based on the evaluation of 
temperature measured by a thermometer in 
contact with the coolant. Given that the dynamic 
response of the thermometer was not 
compensated, results only evidenced the different 
temperatures obtained with liquid nitrogen and 
Dimethyl Ether therapies . Neither the comparison 
of effects on patients nor the measurement of 
temperature of the boiling fluid provide for 
quantitative parameter describing the heat that the 
cooling fluid removes from the skin. In addition, 
both methods prevent from performing a large 
number of tests in repeatability conditions. 

In this paper, we describe a method for the 
measurement of the heat removed by different 
OTC devices for the skin tags removal, which are 
very similar to the ones used for verrucas 
treatment (same fluid, different spray nozzle). The 
measurement method, based on a heat flux meter, 
allows evaluating the heat flux produced by the 
OTC devices. As later explained, the heat flux can 
be used to estimate the skin temperature using the 
skin thermal resistance models available in the 
literature [5,6]. 

2. Proposed Method 

The thermal scheme of the interaction between 
the evaporating gas and the skin is shown in Fig. 
1 a: the heat flux imposed by the medical device 
(Jgas) is partially absorbed by the body (J S ki n ) and 
partially from the environment ( Ji ea k) • Jieak is 
generally small and depends on the OTC device 
characteristics. The skin temperature - one of the 
crucial parameters for treatment effectiveness - 
depends on the gas boiling temperature, on the 
skin thermal capacity and resistance and on the 
leak resistance. The comparison between different 
OTC devices can be performed by creating a skin 
simulator, i.e. a device that has a thermal 
resistance and a thermal capacity similar to that of 
the human skin. In this case, the comparison 
could be based on the cooled skin temperature or 
on any other parameter. An alternative approach 
(the one adopted in this paper) consists in the 
identification of the heat flux actually produced 
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by the boiling gas and the estimation of the skin 
temperature using one of the models existing in 
literature. 

With the proposed measurement method, shown 
in Fig. 1 b, the human body is substituted by the 
heat flux meter, which derives the heat flux by 
dividing the temperature difference measured 
across a thermal resistance (A7fm) by the 
resistance itself (Rfm)- 



Fluid BoilingTemperature 



Interface Temperature 



Body Temperature 



Environment 



Fluid BoilingTemperature 




E ) 



O) 



Environment 



0 



Fig. 1 Thermal schemes of the interaction 
between the human body and the medical 
device for cryotherapy (a) and between the 

device and the measurement system (b) 

Also in this case, part of the flux generated by 
the evaporating gas is absorbed by the 
environment (Ji ea k)» but there are three main 
differences from Fig. 1 a: 

• part of the heat flux entering the instrument 
leaks towards the environment without being 
measured; 

• the heat flux, before being measured, crosses 
a capacitor (C FM ) that modifies determines 
transient response of the measurement device; 
and 

• a thermo electric cooler (TEC) has to be used 
in order to stabilize the temperature of the 
flux meter. 



The heat flux useful for skin treatment purposes 
is Jgas-Jieak- In the scheme of Fig. 1 b this quantity 
equals JpM+JieakFiux+JcFM- JieakFiux has been 
numerically evaluated using the measurement 
model presented in ref. [5] and is shown in Fig. 2. 




Interface Temperature [°C] 



Fig. 2 Heat flux dissipated by radiation and 
convection not measured by the flux-meter 

The lowest interface temperature reached in 
our tests was 5°C. Since the measured J FM was 6 
W, the flux not measured by the flowmeter was 
approximately 500 times lower than the measured 
one (Jieakfiux equal to 0.03 W). 

2.1. Measurement Uncertainty 

According to the model in Fig. 1 b the difference 
between J gas and Ji ea k equals the sum of Jfm and 
Jieakfiux (C F m only affects the dynamic response but 
the net flux on it is zero). Jfm has been computed 
as the ratio between AT FM and R FM . The latter was 
measured using a differential thermocouple, 
calibrated versus Class 1 PT 100 in a thermal bath. 



: AT FM /R FM + Ji e , 



.(Tint) 



(J) 



The measurement uncertainty has been estimated 
by propagating the experimentally identified 
uncertainties in equation 1 using Monte Carlo 
simulations. Given that Jieakfiux is negligible, the 
assumption of an adiabatic model does not 
increase the measurement uncertainty. 




InterfaceTemperature [°C] 



Fig. 3 Relative measurement uncertainty of 
the adiabatic model (squares) and of the full 
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model (circles). Both curves are plotted with 
confidence interval 95%. 

Results evidenced that the measurement 
uncertainty is always lower than 10%. During our 
tests, the expected interface temperature is 
between 10 and 18 °C; consequently, the typical 
uncertainty values range between 2 and 6 % 
(confidence interval 95%). 

Thermocouples were individually calibrated 
using a thermal bath; their standard measurement 
uncertainty was lower than 0.2 °C in the whole 
measurement range, as evidenced in ref . [5] . 

2.2. Dynamic response 

As evidenced in the scheme of Fig. 1 b, the heat 
flux meter has two thermal capacitances, Cfm and 
C H s- The two capacitances reflect in two different 
time constants x, whose values also depend on the 
resistances R FM and R H s- The calculated flux- 
meter capacity is 0.19 J/K. The expected flux- 
meter time constant is 1.25 s and is dominant 
during the cooling phase; given that the typical 
application time is in the order of 30 s, the value 
is satisfying. The second time constant is 
dominant during the heating phase and, as later 
explained, is different for devices with different 
heat flux time histories. 

Although it is theoretically possible to 
compensate the flux-meter time constants (with 
the approach described in references [7, 8]), all 
the tests of interest were performed in quasi-static 
conditions, given that the cooling was much 
larger than the estimated time constant. 

2.3. Tests description 

Nineteen over-the-counter devices have been 
compared using J FM and AT FM as figures of merit 
for the capacity of removing the heat from the 
skin. Three types of device were tested. With the 
"foam" type, the liquid dimethyl ether and 
propane mixture is released into a foam pad 
attached to the device. This foam pad is then held 
on the wart for approximately 20 s. With the 
"spray" devices, dimethyl ether and propane are 
pressurized in a can and are sprayed on the skin 
continuously. With the "dispenser" types, a fixed 
quantity of dimethyl ether and propane aerosol is 
spayed on the skin at each push. 




Fig. 4 Some of the devices that underwent our 
tests. Dispenser type (first, third and fourth 
from the left), foam type (foam not shown, fifth 
from the left) and spray (second and last). 



Aim of the tests was the comparison of different 
nozzles for the dispenser devices. The nozzles 
differed mainly for the airflow passage section 
and location, which affect the evaporation time of 
the dimethyl ether and propane. Placing the 
evaporation holes close to the skin surface 
increases the cooling speed but may increase the 
fluid loss. Conversely, when the evaporation 
holes are far from the fluid, the cooling process is 
slower but the fluid leak is expected to be 
negligible. The different nozzles used in our tests 
are shown in Fig. 5. 



■ - - ' V- 



Fig. 5 Different nozzle specimens. 

For each measurement on a medical device, 
three repetitions were performed. Before each 
test, the interface temperature was restored to 22 
°C using a thermo electric cooler (TEC, Fig. 1 b) 
with the cold side in contact with the heat sink. 
This prevented the progressive cooling of the 
instrument. The temperature data were acquired 
by an Agilent 34970A multiplexer with a scan 
rate of 1 sample/s for each channel. Data were 
acquired via serial interface by a fit-to-purpose 
Lab VIEW virtual instrument. The heat flux and 
the interface temperature were stored on a PC and 
offline analysed. 



3. Results 

The T^x and J FM time histories of the medical 
device number 11 are shown in Fig. 6 and Fig. 7. 
The plots show that the cooling time is 
approximately 10 s, i.e. much larger than the flux- 
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meter time constant (1.25 s). The heating time, 
governed by the second time constant, is 
approximately 40 s. Since the three test are 
performed by actuating the device three times, the 
test repeatability is limited, as evidenced by the 
measurement standard deviation summarized in 
Table 1. 



Table 1 Summary of the tests results 




Fig. 6 Time history of the interface 
temperature reached during the test of the 
device 11_D1 (three repetitions). 




Fig. 7 Time histories of the cooling power 
generated by the device 11_D1 (three 
repetitions) 

The results of all the tests are summarized in 
Table 1. Results allow comparing the different 
solutions in terms of cooling capabilities as 
maximum flux and flux integral. 
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Results show that the spray devices ensure the 
lowest interface temperatures (approximately 17 
°C). The foam devices, thanks to the long 
evaporation time, have a very large flux integral 
in spite of the low AT^. The "dispenser" devices 
are characterized by a AT FM between 3 and 10 °C. 
Given that devices with IDs from 5 to 19 have the 
same dispenser, the role of the nozzle is crucial. 
Best results were obtained with the dispenser 10 
and 11, with 0.5 mm holes placed at 5 mm from 
the skin surface. 

4. Discussion and conclusions 

In this paper, we have proposed a method for the 
comparison of different over-the-counter devices 
for skin tags treatment. Devices using different 
cooling methods (sprays, foam and dispenser) 
have been compared using a heat-flux meter 
originally designed for the measurement of the 
skin thermal resistance. The instrument 
uncertainty in heat flux measurement is lower 
than 6 % (CI. 95 %) and the interface 
temperature uncertainty is lower than 0.3 °C. 

The heat flux produced by 19 different OTC 
devices have been evaluated in controlled 
conditions. Results allowed the identification of 
the best nozzle geometry, which was the one that 
allowed reaching the largest AT FM (or, 
equivalently, the lowest interface temperature) 
and the largest heat flux. 
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Abstract 

In this paper a methodology is described 
concerning the calibration of three-axis low-cost 
accelerometers in the 0 to 10 Hz frequency range, 
to be used for evaluation of existing civil 
infrastructures. Particular attention is paid to the 
evaluation of the calibration uncertainty with 
reference to both static and dynamic operating 
conditions and to the aspects specifically related 
to the test bench and to the calibration procedure. 
These aspects mainly concern the test bench 
behavior in realizing the requested motion law, 
the choosing of the reference for calibration and 
the data processing techniques. 
The obtained results show that a great 
improvement of the low cost accelerometers' 
metrological characterization could be achieved 
according to the procedure described in this 
paper; furthermore, interesting considerations 
have been carried out with reference to the 
evaluation of main and cross sensitivity of 
accelerometers also in dynamic working 
conditions. 

1. Introduction 

Social resilience to disasters is now considered 
a topic of highest political and technical concern 
in advanced nations . 

Social safety and building heritage along with 
sustainable urban development are important 
issues in setting guide lines which aim to improve 
this topic. For this purpose, evaluating the 
conditions and performance of existing civil 
infrastructures is a crucial aspect that allows the 
decision-makers to configure the best lines of 
activity in order to improve the quality of life [1] . 

To have integrated procedures involving 
geotechnical and structural aspects finalized to 
buildings' diagnostics, a distributed sensor 
network is needed, allowing us to operate in a 
selective manner and in the most critical 
situations. This implies the need of a greater 



number of measuring sensors, for both the three- 
axis vibration as well as inclination 
measurements. 

This innovation, together with the need of 
covering wide areas through this multi-sensors 
network in order to include different buildings, 
sets the requirement of lower-cost solutions that 
are nevertheless able to ensure the requested 
uncertainty of measurements . 

Most of the above mentioned measurement 
and cost requirements could be satisfied by the 
well known micro-electro-mechanical systems 
(MEMS) technology and some examples of its use 
for these applications could be found in literature 
[2], even though careful attention should be paid 
to many possible causes of errors [3], especially 
when low-cost and low-precision MEMS 
accelerometers are considered. Due to these facts, 
many studies can be found in literature, referring 
to calibration of these sensors, using mechanical 
reference quantities [3], [4], or using fully 
electrical methods to estimate the sensitivity of 
capacitive MEMS accelerometers in batch 
fabrication [5] . 

Even if different aspects are studied, like bias 
correction and main and cross sensitivity 
evaluation, an exhaustive description of 
uncertainty causes can not be found, in particular 
with reference to the dynamic behavior in the low 
frequency range (0 to 10 Hz). 

Market and literature solutions [2] , [6] , [7] and 
characteristics of the calibration test bench 
designed for this specific application [8], [9], [10], 
[11], [12] have been considered, but the provided 
solutions do not appear completely satisfactory if 
the right trade-off between uncertainty of 
calibration on one hand, and the cost of the test 
bench and the calibration procedure duration on 
the other hand are considered. 

Furthermore, bearing in mind these 
considerations, the first aspect to be considered is 
the limit that could be achieved with reference to 
the uncertainty of sensors to be used for 
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acceleration and inclination, taking into account 
cost requirements as well. 

This goal mainly involves sensor behavior, its 
power supply and conditioning, data acquisition 
technique, and initial and periodical calibration, 
especially in the low frequency range of 
vibrations (0 - 10 Hz). 

This paper aims to investigate the main 
technical and procedural aspects to be considered 
for the definition of a calibration procedure 
modulated in function of the application, as well 
as the cost range of the sensors. 

Particular attention is paid to the experiment 
design, the way through a mechanical reference 
acceleration is created and evaluated and how data 
processing techniques influence the variability of 
results, retrieving useful information without 
having to perform superfluous tests, according to 
the quality level of the sensors that are being 
selected. 

It is to be pointed out that straight through the 
paper the acceleration will be measured in terms 
of g, acceleration of gravity; the authors are aware 
that is formally incorrect, but they prefer this 
solution for better understanding of results in the 
field of civil and seismic engineering. 

2. Methodology 

Scientific evaluations, experts' advice and 
market analysis led us to define the following 
ranges of interest for tests: 

• frequency range: 0 - 40 Hz; 

• maximum amplitude of vibration: ±2g. 

In this paper the calibration in the frequency 
range 0 to 10 Hz will be studied. It is external to 
the best vibration frequency range of standard 
electro-dynamic shakers. 

Calibration uncertainty of sensors will be 
evaluated with reference to the contribution of the 
effects tied to both the calibration test rig and the 
sensor characteristics. 

The three-axis sensor will be calibrated both 
statically and dynamically. In fact, in field 
applications we are interested to use it as an 
inclinometer, the steady behavior, and as a three- 
axis accelerometer, involving the dynamic 
behavior. 

For the general calibration of the 
accelerometer the following relationship between 
the input acceleration and the output signals can 
be set: 



a x \ (S xx S xy S xz \/v x \ /q x \ 

a y) = [ S yx S yy S yz \\Vy\ + \q y \ Eq. 1 

a z J \S ZX S zy s zz )\v z ) W 

being S xx> S yy> S zz the main axis sensitivities; the 
other elements of the matrix are the transverse 

(cross) sensitivities. I a y I are the measured 
acceleration components, I V y I are the voltage 

W 
( qx \ 

outputs of the sensor, and I Qy I is the offset 
\R Z J 

vector. 

A total of 12 parameters have to be evaluated 
[3], [13], [14], [15]. 

For the static calibration a mechanical arm is 
used in order to adjust the sensor to different 
positions, through the gravity acceleration vector 
components along the three axes. 

The minimum number of measurement 
orientations is 4, although more measurements 
may give a more robust calibration, using a linear 
Least Squares Optimization. In this article, the 
effect of the choice of the number of the angular 
positioning on the calibration uncertainty is 
investigated in order to provide useful information 
about the calibration procedure that is more 
suitable for the characteristics of the selected 
accelerometers. 

The behavior of the main and transverse 
sensitivities in the frequency range 0 to 10 Hz is 
also a point of interest. 

For the dynamic calibration a test bench based 
on a rotary device is used. It is driven by a 
brushless servomotor, controlled by a PLC by 
means of a high accuracy angular encoder, which 
allows us to realize different motion laws 
(sinusoidal, saw-tooth, ramp, etc.). 

Being the rotation axis vertically placed, all the 
three axes will operate, measuring the acceleration 
of gravity g (or a component of it, depending on 
the position), the tangential acceleration, a t (or a 
component of it, depending on the position), and 
the centrifugal acceleration, a c , respectively. 
These acceleration values will be used as the 
reference accelerations [5], according to the 
following equations, where CO is the angular 
velocity and r is the radius of sensor positioning: 

a t = (b ■ r Eq. 2 

a c = a) 2 ■ r Eq. 3 
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It is to be pointed out that a t is obtained 
through twice differentiation of angular encoder 
signal, which could cause high frequency noise 
problems. A high accuracy three-axis 
accelerometer will also be employed for further 
validation purposes, according also to literature 
indications [12]. 

The test bench is depicted in Fig. 1 . 

Furthermore, the accelerometer is placed so 
that X-Y or X-Z planes do not correspond to the 
horizontal one, as Fig. 2 shows: in this way all the 
measuring axes will be subjected to a variable 
acceleration at a frequency depending on the 
repetition rate of oscillations. 




Fig. 1: The rotary table, driven by the brushless 
motor. 




Fig. 2: Scheme of tangential acceleration 
components corresponding to a rotation angle a of 
the plane X-Z with respect to the horizontal one. 



In order to evaluate the effect of this parameter 
and to build the dynamic calibration matrix, the 
orientation of the accelerometer under 
investigation will be changed. 



As to the uncertainty of radius, tests are carried 
out at different radial positioning of the sensors, to 
the purpose to reduce the whole uncertainty on the 
reference values of acceleration a t and a c . The 
radial positioning has been varied, inserting 
blocks of increasing and accurately measured 
thickness; the effective position of the rotation 
axis is estimated by a best fitting of the found 
tangential and radial accelerations at different 
radii, hypothesizing a linear behavior with the 
radius. 

The effect of the actual motion law will also be 
analyzed, in relation to its theoretical setting, as 
well as the rotary device set-up. In particular the 
behavior of the calibration bench will be 
examined in case of both square wave and 
sinusoidal motion law for periodical oscillations . 

As for uncertainty evaluation, the following 
assumptions have been made: 

where o r stands for the linearity standard 
deviation from the best fit line (a rectangular 
distribution has been assumed); r = 95.0 mm, o r = 
0.1 mm; 

The contributions to the whole uncertainty of 
a c and of a t deriving from the co and of the (b 
terms respectively, will be evaluated, taking into 
account the repeatability standard deviation of 
them, in correspondence of each angular position. 
Systematic effects on the angular position have 
been considered negligible. 

For geometrical angles a standard deviation of 
1° has been assumed. 

Tests have been planned in order to consider 
the following aspects with reference to the 
uncertainty of a c and of a t (or a component of 
them depending on the position): 

• motion laws; 

• rotary device setting; 

• choosing of the reference; 

• number of points for calibration; 

• data proce s sing technique s . 

Finally it is to be pointed out that the requested 
calibration procedure and the accuracy to be 
achieved should be able to fit the sensor 
characteristics, so that the performances of 
calibration rig and cost of calibration can be 
adequately reduced. 
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3. Results 

3.1. Static calibration 

As a first step, a static calibration of a three- 
axis accelerometer has been done, according to 
the procedure based on Eq. 1 . 

A total of 10 different angular positions have 
been realized in order to fulfil the accuracy 
requirements described in [14]. 

The results are described in the following 
calibration matrix: 



-0.994 
0.0199 
-0.0713 



0.000620 
-0.987 
-0.00310 



0.00404X q 
0.00347 [-] 
-0.994/ 



and the offset vector: 



2.48\ 
2.45 [g] 
2.72/ 



A maximum relative deviation from the mean 
values of 2% has been evaluated for the 
coefficients of the main sensitivities, taking into 
account the angle uncertainty. 

Anyway, measurements have been carried out 
after calibration in the range ±g (±n/2 as for the 
angle) for each axis: a maximum deviation of 
0.0 lg has been evaluated. This result suggests that 
the sensitivity matrix should be evaluated on the 
whole, due to the compensation effect of the 
variability of the different matrix coefficients 
between each other. 

3.1. Dynamic calibration 

3.1 .1 . Identification of the motion 
law 

In case of dynamic calibration of 
accelerometers, as a preliminary matter, the effect 
of the motion time law for oscillations on results 
accuracy has been investigated: in particular saw 
tooth and sinusoidal velocity profile of 
oscillations were considered. 

At a preliminary evaluation saw tooth velocity 
profile, that means square wave profile for a t , 
seemed interesting, having a plateau in a t , being 
very simple to realize on the test bench and 
making able to investigate in the same test many 
frequencies (the odd harmonic components). 

This motion law has been discarded, since 
experimental tests indicated that it is difficult to 
be realized not only practically but also 
theoretically, being necessary specific settings, 
which are meaningless from a physical point of 
view. 



Furthermore, the graphs of Fig. 3 shows that 
strong oscillations in the test bench arise when a t 
is changed, which is a problem also for sensors to 
be calibrated. 




2 3 
Time [s] 

- real tangential acceleration theoretical tangential acceleration 



Fig. 3: Comparison between real and theoretical 
tangential acceleration (saw-tooth motion law, 
fundamental component = 1.6 Hz). 

All the results hereinafter refer to tests driven 
with a sinusoidal motion law; the set-up of the 
rotary table for sinusoidal oscillations is 
remarkably more difficult from a programming 
point of view, since an electronic cam has to be 
realized. Tests have been carried out at a 
oscillation frequency of 0.5, 1.0 and 1.6 Hz; in 
particular indication will be stressed referring to 
1.6 Hz oscillation, whose vibrations peaks are the 
highest. 

3.1.2 Identification of the 
reference 

In a calibration procedure the definition of the 
reference is a fundamental issue. 

The graphs of Fig. 4 show a comparison 
between real and theoretical time behavior of the 
acceleration components, with reference to the 
axes of the accelerometer. It is to be pointed out 
that the components of a t act along the X and Y 
axes, while a c acts along the Z one. 



— real reference x 
— real reference y 
— real reference z 
— theoretical reference x 
— theoretical reference y 
— theoretical reference z 




Time [s| 



Fig. 4: Comparison between real and theoretical 
time behavior of the acceleration components, with 
reference to the axes of the accelerometer. 
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The mean square value of the difference 
between real a^t), a yr (t), a zr (t) and theoretical 
values a xt (t), a yt (t), a zt (t), during the oscillation 
period, has been evaluated, being 0.029, 0.051, 
0.012 g respectively. 

The repeatability of a xr (t), a yr (t), a zr (t) have 
been calculated during the oscillation period. The 
mean values are 0.051, 0.089, 0.016 g, 
respectively. The maximum values are 0.084, 
0.15, 0.030 g, respectively. 

Obviously, most of the oscillations found in 
a xr (t) and a yr (t) depends on numerical 
differentiation, even though higher order 
harmonics can be identified with respect to the 
theoretical behavior, due to effective oscillations 
at points where the motion is inverted. 

Anyway, the relevance of these aspects is 
strongly reduced if the calibration of a three-axis 
accelerometer is considered. 

3.1.3 Data processing techniques 

In Tabb. 1 and 2 D x , D y , D z , are evaluated, 
representing the mean square values of the 
differences between reference acceleration and 
corrected sensor data during the oscillation period. 
In Tab.l reference data are those obtained by 
twice differentiation of angular encoder. In Tab. 2 
the theoretical sinusoidal law is used as the 
reference. The effect of building calibration 
matrix with an increasing number of points, n, is 
also studied. 



n 


D x 




D z 


14 


0.0877 


0.114 


0.0333 


21 


0.0883 


0.112 


0.0299 


28 


0.0829 


0.109 


0.0294 


35 


0.0790 


0.105 


0.0293 


42 


0.082 


0.110 


0.0298 



Tab. 1: D x , D y , D z using the dynamic calibration 
matrix as obtained from the real data. 



n 


D x 


D v 


D z 


14 


0.0780 


0.110 


0.0452 


21 


0.0720 


0.090 


0.0360 


28 


0.0724 


0.0940 


0.0334 


35 


0.0707 


0.0953 


0.0327 


42 


0.0729 


0.0966 


0.0328 



Tab. 2: D x , D y , D z using the dynamic calibration 
matrix as obtained from the theoretical data. 

The differences between the use of the 
dynamic calibration matrix as obtained from the 
real data with respect to the theoretical one are 
negligible. 



Furthermore, provided that a sufficient number 
of points is considered, no effect of increasing it 
has been acknowledged. 

On the contrary, the improvement of using 
dynamic calibration matrix with respect to the 
static one is clearly shown in diagrams of Figg. 
5. a, 5.b, 5.c, where real encoder data, twice 
differentiated, used as the reference, sensor data 
corrected by static calibration matrix and data 
corrected by dynamic calibration matrix are 
compared. 



Fig 5 .a 



— real reference x 

— dynamically corrected 

acceleration x 
— statically corrected 

acceleration x 




Fig5.b 



— dynamically corrected 

acceleration y 
—statically corrected 

acceleration y 
— real reference y 



Fig5.c 



— dynamically corrected 

acceleration z 
—statically corrected 

acceleration z 
— real reference 

acceleration z 



Time [s] 



Fig. 5: Comparison among real encoder data, twice 
differentiated, used as the reference, sensor data 
corrected by static calibration matrix and data 
corrected by dynamic calibration matrix. 

a. a t component (accelerometer X axis) 

b. a t component(accelerometer Y axis) 

c. a c acceleration (accelerometer Z axis) 
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3.1.4 System set up and number of 
points for calibration 

The tests so far carried out also suggested 
some actions to be carried on to improve the 
calibration procedure performances by reducing 
calibration uncertainty, in particular of a t and its 
components; this objective will be pursued by 
making more accurate and repeatable the real 
acceleration profile. 

These actions on one hand require a 
mechanical improvement of the rotary table, 
allowing us to increase the transferability of the 
angular encoder information to the position of the 
accelerometer to be calibrated, on the other hand 
ask for a higher sampling rate of the angular 
encoder data. 

It is to be pointed out that the used PLC allows 
us to use a task time of less than 1 ms, allowing a 
sampling rate of about 1 kHz. 

Rather than increasing the number of points on 
which the matrix should be evaluated but 
provided that a sufficient number of points is 
considered, it seems more convenient to mix data 
measured changing the sensor orientation, also in 
dynamic calibration. 

4. Conclusions 

In this paper some aspects concerning the 
calibration uncertainty of three-axis low-cost 
accelerometers for diagnostics of civil buildings 
are considered. Both static conditions and 
dynamic behavior in the range 0 to 10 Hz as for 
vibrations have been examined. 

A frequency varying three-axis field of 
acceleration has been realized by modulating both 
gravity acceleration, and tangential/radial 
acceleration components of a rotary motion 
created on a test bench driven by a brushless 
servomotor. 

Reference values for a c and a t acceleration (or a 
component of it depending on position) have been 
deduced from the indications of a high precision 
rotary encoder used for axis control, in order to 
separate the effect on the whole uncertainty of 
calibration of the test bench characteristics and 
settings from the sensor performances . 

The effects taken into account refer to : 

• motion laws; 

• rotary device setting; 

• choosing of the reference; 

• number of points for calibration; 

• data processing techniques . 
Experimental results, for both static and 

dynamic applications, show that both main and 
cross sensitivities of a three-axis capacitive 
accelerometer of good quality level from a 



metrological point of view can be obtained; this 
allows to correct the effect of cross sensitivities 
also in dynamic applications, in order to 
remarkably improve the transducer accuracy for 
in field applications . 

The work carried out it is to be considered as a 
first step of the whole research; in fact the 
experimental results have also suggested actions 
to be realized in order to improve the mechanical 
behavior of the test bench, the acquisition of the 
angular encoder data, the calibration procedure 
and the data processing techniques. Improving 
these aspects is expected to further reduce the 
calibration uncertainty. 
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Abstract 

The measurement of contact pressure of human 
fingers is very important to understand human 
perceptual mechanisms, that is the main goal of 
most of the neuro scientific studies. It may also 
lead to a correct development of tactile devices 
and haptic systems, as they are intended to convey 
controllable and effective stimuli. 
In this work, an optical measurement system 
based on Frustrated Total Internal Reflection 
(FTIR) is proposed for the measurement of the 
pressure distribution on the contact area between 
a human finger and a flat surface. The feasibility 
study performed shows that the tested sensor can 
be effectively used for the measurement of the 
fingertip contact pressure both on static and 
dynamic conditions. 

Keywords: Frustrated Total Internal Reflection 
(FTIR) , fingertip contact pressure, biomechanics 

1. Introduction 

The investigation of the mechanisms of human 
tactile perception represents a fundamental topic 
in haptics. The fingertip deformation is the basic 
mechanical action in the tactile perception, since 
tactile sensitivity depends on the tissue strain and 
hence on the contact area. Different models are 
available to predict the behavior of the fingertip in 
terms of contact area, deformation and pressure 
distribution [1,2]. Measurements are thus required 
both to determine the static fingertip area and 
force as input to the numerical models and to 
measure the area and the pressure necessary for 
numerical model validation. When measuring the 
contact pressure, the goal is to set up a 
measurement system with high sensitivity at low 
pressure, since, for specific applications, the 
typical pressure range is from 0 Pa to 50000 Pa. 

The main objective of this work is the 
development of a measurement system, which 
allows to characterize the different aspects of the 
interaction finger - contact surface, through the 



measurement of the contact area and of the 
pressure distribution, both in static and dynamic 
conditions. 

Many sensors based on different physical 
principles have been proposed for the 
measurement of pressure at contact surface 
between two rigid or flexible bodies. Some 
contact pressure measurement systems have been 
used for measuring pressure distribution at 
interfaces between objects, with application to 
comfort analysis and improvement of vehicle 
seats [3] , for various biomedical applications such 
as measurement of contact pressure between hand 
and handle of a tool in order to analyze vibration 
transmission to the hand [4,5] or between foot and 
ground for plantar pressure analysis [6,7]. Many 
applications have been developed [8] also in other 
fields such as robot technology (tactile sensors), 
for contact force mapping of mechanical parts. 
The film sensors used for those applications are 
based on piezo-polymers, capacitive, conductive- 
ink or resistive polymers sensing elements. 
However, many metrological problems remain 
unresolved, due to non-linearity, reological 
behavior, mounting surface curvature and shape, 
dynamic characteristics, etc. Also frequency 
limitations can be due to the sensor itself or to the 
electronic processing and data acquisition system, 
when a matrix of many sensing elements is 
necessary. Piezo-polymers have probably the best 
metrological characteristic but cannot measure 
very slow fluctuations or static pressure 
components. Sensors based on conductive ink or 
resistive polymers showed from preliminary tests 
that metrological characteristics change with time 
and have large hysteresis. The capacitive sensing 
elements can be used to measure both dynamic 
and static pressure with acceptable linearity and 
hysteresis, but their sensitivity is low and they 
cannot measure very small pressures. 

Here an optical measurement system based on 
light scattering is developed and analyzed. 
The outline of this paper is as follows. First, the 
design of an optimized measurement system 
based on Frustrated Total Internal Reflection 
(FTIR) is described. Then the results of specific 
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tests performed on the measurement system are 
shown, to define its static and dynamic behavior 
and to evaluate if the proposed measurement 
technique can be effectively used to measure 
human fingertip mechanical properties in terms of 
contact area and pressure distribution. 

2. Experimental background and 
setup 

The principle of the measurement system is the 
light reflection/refraction at the interference 
between two transparent media. When the light 
passes through an interference between media of 
different refractive indices, the path of travel 
changes (Snell's Law). 

Thus when a light source is installed at the edge 
of a thin glass plate and the incident angle $i 
exceeds the critical angle 0 critica i (Fig. la), the plate 
acts as a light guide. In this case (total internal 
reflection) the light cannot be seen when 
observing the plate from below the plate. 




(a) 



between 'soft' structures and a rigid planar 
surface, specifically for application to tire 
footprint analysis [9,10]. These systems are 
optimized to measure tire footprint pressures and 
are not suitable for the measurement of low 
pressures. The design of the proposed 
measurement system was developed to permit the 
measurement of low pressures that characterize 
human fingertip contact. The roughness and the 
mechanical characteristics of the contact surface 
are the most critical parameters for the realization 
of an effective pressure measurement system 
based on FTIR. These parameters must be 
carefully evaluated, because of their effects on the 
range and on the sensitivity of the measurement 
system. 

The measurement system setup is shown in 
Figure 3: a sliding system realized through 
threaded rods permits to adjust the camera 
position depending on its optics; a 200 mm x 200 
mm x 200 mm darkroom was realized to avoid the 
interference of light from external environment; a 
150 mm x 150 mm x 10 mm glass plate was fixed 
on a square support with a green LED strip 
installed around the perimeter. 

The main characteristics of the glass plate are 
shown in Table 1 . 




Fig. 1: Principle of measurement [9] 

If a reflecting media is pressed against the plate, 
the surface of the media will scatter light in the 
contact surface, being the media with a refractive 
index different than the refractive index of the air 
(Fig. lb). By observing the light scattered, it is 
thus possible to view the shape of the contact 
surface. Depending on the roughness of the 
contact surface, the contact area will increase with 
the contact pressure, thus the light intensity can be 
used to measure the contact pressure between the 
two media (Fig. lc). 

Methods based on white light refraction are used 
for the measurement of normal contact stress 



Chemical glass composition 


(Si0 2 ) + oxides 


Aspect 


Transparent 


Christallinity 


Amorphous 


Refractive index 


1,458 - 1,86 


Nominal tensile strength 


4MPa 


Nominal compression strength 


1 GPa 


Nominal bending strength 


40MPa 



Table 1: Glass plate characteristics 



3. System characterization 

3.1. Static characterization 

A first static analysis of the system was 
performed to characterize the behavior of the 
sensor when a constant pressure is applied on it. 
For this purpose, the glass plate was covered with 
a thin membrane and fixed on a square support 
with a green LED strip installed around the 
perimeter. 

A Plexiglas frame was then realized to install the 
thin lattice membrane at the top of the glass plate: 
the membrane is stretched, thus there is a thin air 
gap (no contact) between the membrane and the 
glass plate in the starting configuration. The 
lattice membrane is talc -covered to increase the 
roughness of the contact surface and thus the 
sensitivity of the measurement system. 
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Figure 2 illustrates a typical result, from the 
literature, of the load-deformation curve for the 
latex membrane used. 

V 

Pressure [N/cm2] 



Thikness reduction [urn] 

Figure 2: Typical pressure - sensor thikness 
reduction curve [11] 

From this result it is possible to see, in general, 
the non-linear behavior of the latex membrane. 
However in the pressure range of interest (0 - 
50000 Pa), as a first approximation, we can 
consider linear the curve pressure-thickness 
reduction. 

For the static calibration, a reference pressure is 
applied to the sensor by the thin membrane, which 
is loaded by a constant pressure chamber, whose 
pressure is measured by a pressure sensor 
connected to a Scandura Pascal 100 system. In 
order to avoid pressure losses, the contact 
between the pressure chamber and the glass plate 
is realized through an O-Ring for pneumatic 
applications (Fig. 3). 



pressure chamber 



inlet valve 
/ 



chamber volume 



v V 

\ O-Ring ' 




O-Ring 




1 / 1 
h n 


lol 




tested area 


glass plate 


0 



Fig. 3: Test bed for static calibration 

The reference pressure varies from 0 Pa to 
100000 Pa, that is the maximum pressure to avoid 
fracture of the glass plate. 

A high performance reflex digital camera 
(Canon EOS HOOD, with a 12,2 megapixel 
CMOS sensor), was used to take high resolution 
grey -level images . The camera was equipped with 
a 18-55mm f/3.5-5.6 lens: the distance from the 
camera to the glass plate was adjusted depending 
on the focal length of the camera. The camera 
parameters (ISO settings, shutter speed, white 
balance) were adjusted to obtain maximum 
sensitivity for the specific setup (f-number = 5,6, 
exposure time = 1/8 s, ISO 400). 



The raw images (4272 x 2848 pixels) were 
acquired with the software Canon Eos Utility, 
converted in .tiff format and thus post-processed 
with Matlab toolbox. Images were acquired for 
each pressure step (5000 Pa - 21 steps) in the 
calibration range. 

For a given pressure, the grey level of each pixel 
varies depending on the pixel location and on the 
local characteristics of the contact surface 
(roughness). The uncertainty is less than 1% full 
scale. Figure 4 illustrates a typical static 
calibration diagram: abscissa refers to the 
reference pressure (qi [Pa]); ordinate refers to the 
grey-level measured by the sensor (qo [bar]). The 
calibration curve is obtained using the method of 
least squares. 



Figure 4: Typical static calibration diagram for a 
pixel 

After these first measurements, a measurement 
area was defined to reduce the dispersion due to 
the localized imperfections of the contact surfaces 
(Fig. 5). This is the most effective sensor area 
both for static and dynamic measurements . 




Fig. 5: measurement area selection 
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3.2. Preliminary dynamic 
characterization 

Then a dynamic characterization was performed 
by loading the sensor with a sinusoidal excitation 
applied through a rigid feeler pin (Fig. 6). The 
compression force, generated by the 
electrodynamic exciter, is measured by a load 
cell. 



For the images processing a specific software in 
Lab View has been developed (Fig. 8) 




Fig. 6: test bench for dynamic characterization 

Tests were performed using sine excitation: 
typical signals obtained are shown in Figure 7 . 





Fig. 8: image processing software interface 

Changing excitation frequency keeping fixed the 
level of force the ratio between system output 
(gray level) and force amplitude vs frequency is 
plotted in Figure 9 . 

I nt*nslty/Faix« 

Lgny I* Ye I / Nj 



1? 



Figure 7: Time history output system and force. 



0 b io IS 

Figure 9: Time history output system and force. 

These results allow establishing that a useful 
bandwidth is from 0 up to 10 Hz with tolerable 
response decrease lower than 5%. 

4. Conclusions 

The analysis of the contact pressure between the 
finger and a contact surface is of remarkable 
interest in the field of biomedical design, to 
characterize human fingertip mechanical 
properties in terms of contact area, deformation 
and pressure distribution. 

In this work a measurement technique based on 
FTIR is proposed for the analysis of contact 
pressure between the finger and a flat contact 
surface (glass plate) . 

A FTIR measurement system was designed and 
tested to evaluate the feasibility and effectiveness 
of this measurement technique to measure the 
contact area and the pressure necessary for 
numerical fingertip model validation. 

Results showed that there are good assumptions 
to use the sensor for the measurement of the 
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fingertip contact pressure both on static and 
dynamic conditions. 

Consistency and repeatability observed at this 
first stage is very encouraging for the 
development of an effective fingertip pressure 
measurement system based on FTIR. 

This will permit to analyze the behavior of the 
finger for different conditions, and thus to 
validate numerical models not only in contact area 
and force, but also in pressure, that is suitable for 
investigating the fingertip deformation in static 
and dynamic conditions. 

Plans for future work are to further explore new 
materials for higher sensitivity, test different 
configurations with a complete set of loading 
conditions. 
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Abstract 

An advanced measurement system has been 
developed to provide the real-time monitoring of 
indoor thermal comfort. The core device is a low- 
cost IR sensor assembled with two fixed-step 
motors, installed on the ceiling of the occupied 
room. The embedded tool performs the automatic 
scanning of each indoor surface to evaluate the 
temperature distribution. The mean radiant 
temperature (t r ) and predictive mean vote (PMV) 
are computed for several positions in the 
environment and provided as output of the device 
through wireless or wired connection. The study 
presents the calibration of the IR sensor to allow 
the continuous measurement of the indoor 
surfaces temperature to retrieve the t r . The results 
of the calibration procedure are used to employ 
and validate the IR-based system in a real 
environment, a classroom. The validation entails 
the comparison with a microclimate station and 
the subjective response of the occupants. 

1. Introduction 

According to the rational approach [1], the 
evaluation of thermal environments requires six 
quantities: two personal parameters (clothing 
thermal insulation and metabolic rate) and four 
physical parameters (air temperature, mean 
radiant temperature, air velocity and air 
humidity). These parameters have to be measured 
or estimated in order to derive the indoor thermal 
comfort by means of the PMV (Predicted Mean 
Vote) index [2]. As presented in [3], the mean 
radiant temperature t r is a very significant factor 
in moderate environments, especially in buildings 
whose envelopes are exposed to a strong solar 
radiation [4] . A method based on the computation 
of angle factors, provided by IS07726 [5], is 
capable of estimating the t r with a good accuracy. 
However, it is complex because of the need of 
measuring all the surfaces temperature inside the 
room and the solar load with radiometers in case 
of transparent surfaces. For example, in [6] a 



solution based on an infrared (IR) sensor is 
presented to estimate the thermal comfort of a 
subject in a vehicle in which the mean radiant 
temperature is not measured following the 
standard procedure in [3], but it is assumed to be 
equal to the measured local temperature. This 
paper presents an ad-hoc solution to provide the 
PMV index by means of continuous measurement 
of the indoor surfaces with an IR-based system 
and application of the angle factors methodology 
to calculate the t r . Basics of the approach were 
presented by the authors in [7] with an initial 
validation to demonstrate its feasibility. The final 
system was presented with the description of the 
methodology and validation in a real office in [8] 
and [9]. In [10] a specific application of the 
system in AAL was illustrated. The complete 
description of the calibration procedure of the IR 
sensor was presented in [11]. 
This paper presents shortly the system and the 
calibration of the IR sensor to achieve the 
measuring performance required to allow its 
deployment. Then, the validation in a real 
classroom based on the comparison with a 
microclimate station and the subjective 
investigation about the thermal environment is for 
the first time illustrated to demonstrate the 
reliability of the solution proposed in a practical 
case currently of large interest. In fact, several 
studies concerning the methodology to assess the 
thermal comfort in public buildings are present in 
literature, such as sports and leisure facilities in 
[12], or, specifically for the classroom 
environment, in [13] [14]. The problem of how to 
accurately measure the thermal comfort in 
classrooms is thus a field of investigation where 
the solution proposed can address the issue of 
overcoming the traditional measurements 
performed with a thermostat or a temperature 
probe in the return air duct of the ventilation 
system. 

The measurement system is now patent pending 
[15]. 
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2. Assessment of the thermal comfort 

2.1. Comfort model 

According to ISO 7730 [2] "a human being's 
thermal sensation is mainly related to the thermal 
balance of his or her body as a whole. This 
balance is influenced by physical activity and 
clothing, as well as the environmental parameters: 
air temperature, mean radiant temperature, air 
velocity and air humidity". In a moderate 
environment, the human thermoregulatory system 
will automatically attempt to modify skin 
temperature and sweat secretion to maintain the 
heat balance. As the definition given in ISO 7730 
implies, thermal comfort is a subjective sensation 
and it can be expressed with a mathematical 
model of PMV which is function of four 
environment parameters and two personal ones: 
PMV = f(t a ,t r ,v a ,RHJ cl ,Af ) Eq. 1 

where t a [°C] is the room air temperature, t r [°C] 
is the mean radiant temperature, v a [m/s] is the air 
velocity, RH [%] is the relative humidity, I c i 
[m 2 K/W] is the clothing insulation and M [W/m 2 ] 
is the metabolic rate. The PMV index predicts the 
mean value of the votes of a large group of 
persons with a 7 -point thermal sensation scale 
(from -3, cold, to +3, hot) exposed to the same 
environment. The equation comes from the heat 
balance between the human body and the 
environment. The PPD (Percentage of People 
Dissatisfied) index is used to provide information 
on the thermal discomfort or thermal 
dissatisfaction predicting the percentage of people 
likely to feel too warm or too cool in a given 
environment. The PPD is calculated as function of 
the PMV: 



PPD = 100-95 • e 



-(0.03353^MV 4 +0.2179PMV 2 ) 



Eq.2 



An environment is considered as very 
comfortable when the PMV varies between -0.5 
and +0.5. It is considered as comfortable for 
values between -1.0 and +1.0. These cases lead to 
a PPD of 10% and 27% for -0.5<PMV<+0.5 and - 
1.0<PMV<+1.0 respectively. When PMV is zero, 
that is to say for the perfect case, the PPD is 5%. 

2.2. Standardised measurement 
procedure 

The measurement of comfort indexes is 
generally performed with a microclimate station. 
The station includes the probes needed to acquire 
the ambient parameters and a data logger. Then, a 
software provides the capability to apply the 
personal parameters and calculate the comfort 



index in post-processing. This procedure is used 
for short period monitoring and not for real-time, 
continuous monitoring. The station used for this 
study is the "HD32.1". PMV and PPD indexes are 
calculated with DeltaLog 10 software. The 
metrological characteristics of the probes used are 
recapped in the Table 1 . These characteristics can 
be considered representative of the standard 
measuring performances for standard comfort 
assessment according to the ISO 7726 [5]. 
Table 1: Measurement characteristics of the 
microclimate station 



Range 



Uncertainty 



-30°C 4 
100°C 



-30°C + 
120°C 



t r Globe 

thermometer 

v a Omni directional 0.05 - 
hot wire 

RH Capacity 5 -98 



±0.1°C 



±0.02m/s (0.05-lm/s) 
±0.1m/s (l-5m/s) 



2.3. Subjective investigation 

The subjective approach has the aim of 
investigating the perception of the thermal 
environment by means of judgment, preference 
and acceptability through surveys of the users. 
This was performed using questionnaires, that 
users were called to compile, elaborated from the 
model provided by ISO 10551 [16]. This method 
allows from one side the investigation of the 
occupants perception and from the other side to 
assess the capability of the measuring system to 
evaluate the comfort level and the typical 
constraints to apply a certain measurement 
methodology. An example of the subjective 
investigation is reported in [17]. 

There are a number of subjective judgment 
scales for thermal environments, which differ in 
whether emphasis is placed on some aspects of 
judgment: perceptual or affective (evaluative and 
preferential), global or localized, present or past, 
instantaneous or extended over a period of time. 
The survey drafted and used for this validation is 
shown in Fig. 1. The occupants were asked to 
provide general information (age, sex, kind of 
activity performed and clothing level), mark their 
position in the classroom and then indicate their 
thermal perception (point 6), thermal preference 
(point 7) and affective assessment (point 8). 

The analysis of the thermal perception votes was 
performed according to standard [16]. It asserts 
that the only study devoted to the statistical 
characteristic of data distributions concerns those 
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obtained by applying the 7-degree scale, from -3 
to +3. Degrees other than the two extreme ones 
are physiologically located at equivalent distances 
and may thus be treated as continuous data whose 
differences are defined numerically (interval 
scale); they approximate normal distributions, so 
it is quite legitimate to calculate the mean and 
standard deviation. 
The degrees on the scales of perception proposed 
(thermal preference and affective assessment) are 
ranked as classes of observable data that 
correspond to a totally ordered finite 
mathematical set (ordinal data or scales). 
Therefore, the appropriate statistical tests of the 
null hypothesis are of the nonparametric type (e.g. 
sign test, median test or certain types of variance 
analysis) . 

The analysis presented in this paper focuses on 
the comparison between the thermal sensation and 
the PMV measured. Thus, the thermal preference 
and affective assessment are not considered for 
the comparison. 

Survey on thermal environment - OCCUPANT 



1. Age: 



3. Date: 



4. Level of activity performed: 



5. Clothes within the test (mark the most 
appropriate solution) 

□ Underwear, socks, tracksuit, t-shirt, shoes 

□ Underwear, shirt, trousers, socks, shoes 

□ Underwear, shirt, trousers, smock, socks, shoes 

□ Underwear, shirt, trousers, jacket, socks, shoes 

□ Underwear with long sleeves and legs, shirt, 
trousers, V-neck sweater, jacket, socks, shoes 

□ Other 



Mark your position in the room 



6. How do you feel at the precise moment (mark appropriate box): 





Very 
hot 


Hot 


Warm 


Slightly 


Neither hot 
nor cold 


Slightly 
cool 


Cool 


Cold 


Very cold 


Begin of 




















Break 




















End of 





















7. At this moment, would you prefer to be ... ? 





Much 




Slightly 


Without 


Slightly 


Cooler 


Much 


Begin of 
lesion 
















Break 
















End of 

















8. At this moment, do you find this ...? 





Perfectly 


Slightly difficult 


Fairly difficult to bear 


Unbearable 


Begin of lesson 










Break 










End of lesson 











3. System description and 
methodology 

3.1. Methodology applied 

The solution proposed is designed to provide the 
performances of a microclimate station, but in a 
manner suitable for the continuous real-time 
monitoring and for multiple positions. Fig. 2 
represents the concept of the system. For this 
reason, the system is built with low cost 
components to measure all the needed ambient 
parameters. The capability to measure mean 
radiant temperatures with a spatial discretization 
and with a single sensor represents the main 
advance of the system. The specific design and 
calibration for this scope allows the employment 
of such solution in the field of the thermal 
comfort assessment. 



Central 
Unit 




Fig. 1: Survey used to assess the thermal comfort by 
subjective judgements 



Fig. 2: Concept of the IR scanning system 

The whole procedure to configure, control and 
acquire data is managed by a control unit. It is 
composed of a commercial board (Arduino Mega) 
and all the electronics needed for sensors 
integration and wired/wireless communication 
(i.e. Bluetooth for interaction with Android 
devices through a dedicated GUI). The inputs 
required to perform the algorithms embedded can 
be stored directly into the microcontroller 
EEPROM memory, or updated by the user and/or 
the technician. 

The user can interact with the system through 
the terminal with the user interface implemented 
(PC or Android device) and is asked to provide 
the inputs required for the calculation of thermal 
comfort index (room geometry, subject positions 
and their personal parameters, metabolic activity 
and clothing thermal insulation). The 
microcontroller manages the automatic scanning 
of all the surfaces of the room, provides the 
needed measurements through integrated sensors 
and formulas implemented into the embedded 
C++ libraries and then sends the results back to 
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the device for the real-time monitoring and post- 
processing. 

Next paragraphs describe in details the 
components used and the methodology applied to 
address the specific measurement requirements. 
Further details can be found in [1 1] . 

3.2. The IR scanning system 

The concept relies on an IR scanning system to 
install on the ceiling of the room, that measures 
the indoor surface temperatures, and sends them 
to the control unit to calculate the mean radiant 
temperature t r and PMV for several locations 
inside the room. The mean radiant temperature is 
derived for multiple positions of the subject inside 
the room according to the angle factors algorithm, 
as presented in the ISO 7726 [5]. The t r is 
computed from the weighted average of the 
internal surface temperatures u and the respective 
view factor in relation to a subject, F s -u for N 
surfaces: 

C=2X,7r Eq.3 

i=l 

A mathematical expression developed by 
Cannistraro et al. [18] was used to calculate the 
view factors between a subject and internal 
surfaces. This approach, turned into an algorithm, 
automatically computes the coefficients 
differentiating the kind of surface considered 
(vertical or horizontal) . 

The procedure to measure t r requires the 
measurement of all the indoor surfaces 
temperature. Therefore, the core of the device 
consists of an array of thermopiles for non-contact 
IR temperature measurement and two servos to 
manage the orientation of the sensor (0°-180° in 
both directions), as shown in Fig. 3. This device, 
installed on the ceiling of the room and possibly 
in the center, allows the continuous measurement 
of the indoor surface temperatures. Depending on 
the Field of View (FOV) of the sensor (5.12° by 
6° for each pixels in the row), an embedded 
algorithm was implemented into the 
microcontroller to manipulate the temperatures 
measured iji(h,v)j2(h,v), ...T6(h,v)), where h and v 
represent the horizontal and vertical orientation of 
the sensor during the acquisition. The application 
of this algorithm provides low-resolution thermal 
images as output for each surface. The IR sensor 
is a commercial solution (array of eight 
thermopiles arranged in a row, built in electronics 
and a silicon lens), which provides a temperature 
measurement of eight consecutive points with a 
resolution of 1°C each. A methodology to correct 
the surface temperature, based on its emissivity 



and environmental factors, is applied as part of 
the complete measurement process, as discussed 
above. 



a) ^ b) 




Fig. 3: IR scanning system, a) Assembled device (IR 
sensor + servos); b) Detail of the thermopile array 
adopted 

The issue of correcting the infrared emissivity is 
addressed to reduce the measurement uncertainty. 
In fact, an on-board correction of the IR raw 
measurement is implemented by means of the 
following equation: 

T 4 — _Lt 4 - — 8 T 4 - — T T 4 F A 
1 obj~ 1 tot 1 refl 1 atm 

ST S ST 

where s is the surface average emissivity, z 
is the transmission coefficient of the atmosphere 
(assumed as a constant value of 0.99), T m is the 
total temperature (raw value) measured by the IR 
sensor, T re fi is the reflected temperature and Tatm is 
the temperature of the atmosphere (equals to 
indoor air temperature). The average emissivity of 
the surface has to be provided as input of the 
algorithm, while the reflected temperature is 
computed directly measuring the temperature of 
the opposite surface with an emissivity set to 1. 
Given that the IR sensor is installed on the 
ceiling, its inclination angle could affect the 
accuracy of the measuring procedure. However, 
the maximum angle of incidence of the IR ray is 
generally equal to 25°: according to the polar 
diagrams of the infrared emissivity of the material 
in consideration, as in [19], the emissivity can be 
assumed to be constant within certain limit values 
of the incidence angle. In any case, this effect was 
considered and assessed during the calibration 
procedure (paragraph 4) . 

3.3. Environmental sensors 

The evaluation of PMV requires the 
measurement of air temperature, relative humidity 
and air velocity according to ISO 7726; for this 
reason, commercial low-cost solutions are 
adopted for the scope. 

An integrated T/RH sensor allows the single- 
point measurement of the air temperature (t a ) and 
relative humidity (RH) parameters. It is based on 
DS18B20 1-Wire digital thermometer to measure 
the temperature and polymer humidity capacitor 
for RH and they are both connected to an 8 -bit 
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single chip computer. The output is a calibrated 
digital signal that derives from a previous 
calibration and coefficients saved in the OTP 
memory. 

A low-cost flow sensor (1ST FS5), based on the 
thermal conductometric principle, has been 
adopted for single -point air velocity 
measurement. The sensor includes two platinum 
resistors in a chip and can be used as a constant 
temperature anemometer (CTA) through an 
electric circuit provided by the sensor 
manufacturer. An ad-hoc PCB was developed and 
integrated in the control unit to derive the air 
velocity, as a function of the bridge voltage 
measured by means of a 12-bit A/D Converter. 
The calibration curve has been found by means of 
a calibration conducted in the UNIVPM 
laboratory, with a reference anemometer (probe 
AP3203 of HD32.1 microclimate station). After 
calibration, the sensor presents an accuracy of 
±0.06 m/s in a range of 0-=-l m/s. 

4. IR system calibration 

The IR system was calibrated and validated with 
experiments conducted in a test chamber [11]. 
The test room is a climatic chamber with two 
black-painted surfaces, which can be heated or 
cooled. Surface temperatures, together with 
environmental parameters, can be set according to 
the operative ranges. Four k-type thermocouples 
were placed on one of the heated surfaces (Fig. 4 
a) and were used to calibrate the IR system. 




Thermocouple [°C] 
Fig. 4: Test chamber setup (a) and calibration curve 
of the IR system (b) 



Depending on the FOV (Field Of View) of the 
IR sensor and its position within the test 
conducted, a matrix of 5x8 samples was provided 
as final output for each acquisition, where each 
point is the IR temperature measurement in an 
area of 20x24 cm. The calibration was performed 
comparing the values measured by the 
thermocouples with the mean value of four ROIs 
(region of interest). In particular, those four ROIs 
were identified, where the thermocouples were 
placed and the analysis was conducted for the 
average values of the ROIs. A previous analysis 
was conducted to determine the emissivity of the 
surface, which turned out to be 0.89. A test of 7 
hours was conducted when the surface was heated 
up and cooled down, covering a range of 15°C- 
45°C which is typical of indoor environments. 

The calibration (Fig. 4 b) revealed a sensitivity 
of 1.0, an offset of 4.2 °C and an uncertainty of 
±0.9 °C (coverage factor k=2). 

The calibration results were used to assess the 
overall measurement approach to derive the mean 
radiant temperature inside the test room. Two 
surfaces, SI and S2, were heated and cooled 
following predefined profiles in order to simulate 
a real case situation, for example, the solar 
radiation that heats a wall exposed to the exterior. 
During the test, the IR sensor scanned all the six 
surfaces and the other sensors acquired the 
ambient parameters. Therefore, the real working 
operation was simulated in the controlled 
environment. 
Table 2: Measurement of surfaces temperature: 
comparison with the reference 



Parameter 




IR system 


Mean [°C] 


STD [°C] 


T (front surface) (SI) 


-0.1 


±0.8 


T (rear surface) 


0.3 


±0.4 


T (left surface) 


-0.1 


±0.9 


T (right surface) (S2) 


0.1 


±1.8 


T (floor) 


0.1 


±0.6 


T (ceiling) 


-0.5 


±0.4 



The measurement performance is recapped in 
Table 2. Using those temperatures the t r was 
calculated and compared to the one measured at 
the same time with the globe -thermometer of the 
microclimate station. 
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30 



"24 




0 2 4 6 8 

Time [hours] 

Fig. 5: IR system vs globe -thermometer during the 
test. 

The metrological performances turned out to 
provide an average deviation of ±0.5 °C (k=2) on 
the t r with respect to the globe thermometer. This 
result confirms the capability of the solution 
proposed given that the ISO 7726 calls for a 
required accuracy of ±2°C and a desirable 
accuracy of ±0.2°C. 

5. Validation in a classroom 

Given the calibration results, the system has 
been tested in a real case to validate its 
functionality. The validation procedure consisted 
of a comparison with the microclimate station and 
the subjective investigation. The test was 
performed in a classroom at UNIVPM 
(10x9.2x3.4 m) during a regular lesson. The 
students attending the lesson were asked to 
compile an ad-hoc survey to express their thermal 
sensation and preference. The IR system was 
placed on the ceiling of the classroom, nearly in 
the centre and the control unit, with the other 
environmental sensors, positioned in a desk to 
provide the acquisition of the ambient parameters 
and data storing for post-processing. Data were 
compared to that monitored by HD32.1 
microclimate station in two different positions, as 
shown in Fig. 6. Statistical analysis were carried 
out following guidelines provided by [16]. 



HD32.1 
End of lesson 



Control Unit 
Environmental sensors 



IR sensor 
Ceiling mounted 



HD32.1 £ 
Begin of lesson ™ 



Fig. 6: Scheme of the microclimate station HD32.1 
and the proposed solution during the validation in a 
classroom 



5.1. Comparison with the reference 
system 

The comparison of thermal comfort parameters 
retrieved from both the IR system and the 
microclimate station revealed discrepancies of 
±0.5°C for the mean radiant temperature and ±0.1 
for the PMV. Next two figures show the 
comparison between the IR system and the 
microclimate station. 



°L29 
2 

|27 
E 

£25 

c 
ro 

"§23 



-Microclimate station 
• IR system 



14:30 



14:53 



15:16 
Time [hh:mm] 



15:39 



16:01 



Fig. 7: Mean radiant temperature measured by the 
IR system compared with the reference 
measurement 




Fig. 8: PMV measured by the IR system compared 
with the reference measurement 

The capability of the system to provide a spatial 
discretization is also assessed with the two graphs 
above. In fact, the curves compare from 14:30 to 
15:16 the PMV estimated by the IR system in the 
position of the microclimate station at the 
beginning of the lesson (see Fig. 6) while from 
15:16 to 16:01 in the position at the end of the 
lesson. 

5.2. Comparison with subjective 
evaluation 

During the lesson a total number of 40 students 
were analyzed (14 males and 26 females), with an 
average age of 21 ±3 years, a metabolic rate of 1.2 
met (sedentary activity level) and a clothing level 
of 1.2±0.3 clo. Data related to thermal perception 
assessment are shown in Fig. 9. In all the cases 
people perceived a sensation of slight discomfort 
(average PMV of 0.9±1.0): the significance of the 
difference between observed values Y to the 
neutrality level Y =0 was verified by statistical 
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inference (t test for n<60). Results for all the 
cases are summarized in Table 3. 

20 



iBeginning of lesson 
iDuring break 
End of lesson 




-1 0 1 
Thermal Sensation 



Fig. 9: Thermal perception at beginning of the 
lesson, during the break and at the end of the lesson 

Table 3: Statistics of thermal perception analysis 



Parameter 


Begin of the 
lesson 


Break 


End of 
lesson 


Mean 


1.1 


0.6 


0.9 


Standard deviation 
(K=2) 


0.9 


1.1 


0.9 


tn-1 


7.7 


3.4 


6.3 


p-value 


< 0.0005 


0.005 < p < 
0.0005 


< 0.0005 



/7-values were under the admissible error 
threshold (0.05) in all the cases, so the slight 
discomfort perceived was statistically significant. 
The thermal sensation perceived was then 
compared to values provided by the monitoring 
systems as showed in Fig. 10. Confidence 
intervals were computed and compared to validate 
the comfort measurement system with respect to 
the subjective perception. As shown in Fig. 10, 
the subjective assessment revealed a higher 
discomfort level perceived with respect to that 
provided by both the measurement systems, but 
the differences were not statistically significant. 
This demonstrated that comfort measurements are 
comparable with the thermal perception. The IR 
system can provide information about indoor 
thermal comfort that is in agreement with persons 
perceived sensation. Fig. 11 shows the 
comparison in terms of PPD (Percentage of 
People Dissatisfied) index, which was computed 
by both measuring systems. 



1 

-0.6 
0.4 









I I 







Survey IR system HD32.1 
Fig. 10: Comparison of the thermal sensation 
measured with the IR system, microclimate station 
and subjective survey 



Survey I R system HD32.1 
Fig. 11: Comparison of the percentage of 
dissatisfied measured with the IR system, 
microclimate station and subjective survey 

6. Conclusions 

The results show that the system proposed is 
able to provide the real-time measurement of 
surfaces temperature with enough accuracy to be 
employed for the methodology given by IS07726. 
The system provides spatial information of the 
surface to measure the average surface 
temperature with an uncertainty of ±0.9°C (k=2) 
and derive the mean radiant temperature. In fact, 
the validation showed an uncertainty of ±0.5 °C 
(k=2) on the t r and the consequent uncertainty of 
±0.1 (k=2) on the PMV with respect to 
commercial solutions (microclimate station). The 
uncertainty of the t r obtained with the IR system 
is inside the accuracy required by the ISO 7726. 
The application in a classroom, which is a space 
larger than a typical office room, proved the 
capability of the system to provide a spatial 
information. In fact, the comparison with the 
microclimate station showed the same 
correspondence in two different positions. This is 
an important achievement given the relevance of 
the real-time monitoring of thermal comfort in 
schools. Moreover, the comparison with the 
subjective investigation proved the capability of 
the system to assess the thermal sensation of the 
occupants. The measurement system revealed and 
quantified the thermal discomfort perceived 
during the lesson with a degree comparable with 
the subjective sensation. 

Thus, the solution proposed can be employed 
for the real-time and spatial monitoring to 
enhance the control of indoor environments. 
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Abstract 

The present paper demonstrates the applicability 
of a laser-ultrasonics procedure to improve 
performances of train components ultrasonic 
inspection. The method, previously developed by 
the authors, exploits an air-coupled ultrasonic 
probe which detects ultrasonic waves generated 
by a high- energy pulsed laser. As a result the 
measurement chain is completely non-contact 
from the generation up to the detection, this 
letting speed up the inspection time and make the 
set-up more flexible. In this paper the technique 
developed has been applied to train component 
diagnostics, and specifically to wheels and high 
speed train hollow axles with variable section. In 
addition a sensitivity analysis to the surface 
conditions has been performed. 

1. Introduction 

Inspection of train components is a very 
important procedure which has to be periodically 
performed by train operators. Currently periodic 
inspection of axles and wheels is carried out by 
ultrasonic techniques using, for example, phased 
arrays [1] or in particular for hollow axle 
diagnostics a number of rotating contact probes or 
a bore probes [2] . 

Ultrasonic techniques, which are a long 
established method, have the disadvantage of 
requiring the probes to be in contact with the 
object to be investigated, which lengthens the 
inspection time necessary to prepare the object 
and apply the coupling medium. However, the 
typical problems of conventional techniques can 
be overcome by using a hybrid laser-ultrasonic 
system based on the detection of ultrasonic waves 
generated by a high-energy pulsed laser via air- 
coupled ultrasonic transducers [3, 4]. Air-coupled 
ultrasound inspection has already been 
successfully used in many industrial applications 
e.g. NDT on both thick [5] and thin [6] 
composites, NDT on light thin historical vaults 



[7], density measurement of ceramic tiles [8], 
wood detection [9], and thin metallic laminated 
[10,11]. 

The hybrid system proposed is completely non- 
invasive and it makes it possible to: 

- overcome shape and accessibility problems, 

- avoid the application of the coupling medium, 

- shorten the inspection time for large surfaces. 
Laser-ultrasonics has been applied in the 

aeronautical field for damage detection on thick 
composite materials [12] and honeycombs [13] 
and in the railway field for rail [14], rail wheels 
[15] and axle inspection [16]. Gonzales et al. [16] 
presented a Laser Air-Hybrid Ultrasonic 
Technique based on air coupled ultrasounds. The 
limits of the technique are that it is working in 
ablative regime and it exploits Rayleigh waves 
reflected by the crack. The necessity of working 
in ablative regime compromises the surface 
integrity and does not allow considering the 
technique as non-destructive. The use of Rayleigh 
wave reflection makes extremely difficult to work 
with complex geometry, as hollow axles with 
several section variations. In fact, the reflections 
due to the geometry can overlap with the one 
induced by the crack and thus the desired signal 
can be buried into the noise. The present paper 
shows the experimental application of laser 
ultrasonics to the inspection of a real high speed 
train axle and wheels provided by the Italian 
railway company (Trenitalia), where typical 
fatigue defects for the axle [17] and dominant 
wheel failure cracks for the wheels [18] have been 
expressly created according to defects occurring 
in real wheelset [19] . 
For the axles, it has been proved by the authors 
[20] that the technique can be applied also when 
the distance between laser source and ultrasound 
probe exceeds 200 mm. Consequently it is 
possible to apply it for the diagnostics of flaws 
located in the fretting surface by locating the laser 
source and the ultrasound probe at the opposite 
sides of the wheel press fitted area. 
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Fig. 1: Test item (a) and defect positions (b) 



In [21] a model-assisted Probability of Detection 
curve has been derived for this technique, thus 
allowing to define the method reliability. 

Here the attention is focused for the first time on 
the wheels, for which several configurations of 
the laser ultrasonics probes have been tested and 
two significant ones are reported. In the first 
configuration the ultrasound probe is located in 
the thread and works in transmission. In the 
second configuration, the harshest one, the laser 
source and the ultrasound probe are settled on the 
opposite sides of the wheel, thus allowing to 
detect the defect presence in the complete wheel 
profile. 

The paper is structured in three Sections. Section 
2 presents a feasibility study on a train hollow 
axle with section variation. Section 3 describes 
the application to the wheel test case. Finally a 
sensitivity analysis to the surface test item has 
been performed and its results illustrated in 
Section 4. 

2. Experimental results on train 
axle 

The application to a real train axle (Fig. la) 
where typical fatigue defects have been realized 
will be shown and results discussed. In practice 
four convex defects were machined on the 
external surface, see Fig. lb, two of them in the 
wheel fitting area (Dl and D3) and the other two 
in the section transition (D2 and D4). The defect 
morphology is reported in Table 1 . 

In this paper only results corresponding to the 
D4 defect will be reported which is the most 
critical case it being located in the transition after 



the wheel fitting area. For a complete analysis the 
reader can refer to [20] . 

Table 1 Defect morphology. 



Defects 




Size of defect 




H(mm) 


L (mm) 


D Max depth (mm) 


Dl 


31 


1.1 


1 


D2 


28 


1.1 


1 


D3 


28 


1.1 


1 


D4 


28 


1.1 


1 



2.1 Test bench 

The train axle was mounted on a support that 
made it possible to control its rotation, in order to 
be able to scan the axle surface along a 
circumference. The probes were installed on a 
frame where they could move along the axial 
direction. Therefore the complete lateral surface 
of the axle could be inspected by the laser 
ultrasonic system, see Fig. 2. 

The laser ultrasonic system was made up of a 
pulsed laser source, a Nd-Yag IR laser (1064 nm), 
pulses of 12 ns duration and 82 mJ energy, from 
Continuum, and a 1 MHz air-coupled ultrasound 
piezoelectric probe from Ultran Group (model 
NCT210, with 8 mm diameter of active area). 

The ultrasound probe conditioning system was a 
DPR 300 Pulser/Receiver from JSR Ultrasonics. 

The ultrasound signals were amplified with a 
gain level of 69 dB and acquired with a high 
speed Digitizer board NI PXI-5122 (100 MHz 
bandwidth). The laser beam was guided towards 
the axle under test by means of an arm connected 
to the pulsed laser cavity as shown in Fig. 2. 
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Pulsed Laser Cavity 



Pulsed laser output 





Fig. 2: Scheme of the laser-ultrasonic scanning system 



The laser ultrasonics system was installed with 
the laser source and receiving probe far apart, in 
order to point outside the wheel fitting surface. 

A collimated laser beam with a diameter of 
about 8.5 mm was used to keep the ultrasonic 
waves generation within a thermo-elastic regime. 
The air coupled ultrasound probe was inclined, 
with respect to the axle surface normal, of 8 deg 
in accordance to the Snell's Law. This inclination 
was optimized for the reception of surface 
Rayleigh waves, that are the most efficient for 
surface defect identification [20] , in particular in 
case of complex geometry objects. The axle was 
made to rotate by an electric motor and the 
rotation angle measured by an integrated encoder. 
A circumferential scan has been performed along 
an angle of 93 deg with an angular resolution of 
1.5 deg. 




Fig. 3: The laser-ultrasonic scanning system 



2.2 Analysis of results 

For each configuration, a series of ultrasound 
time histories with a duration of 500 pis was 
acquired at every scanning position along the arc 
considered. 



The experimental set-up (i.e. laser and probe 
positions), the axle section tested and the defect 
location (D4) are shown in Fig. 4a. The time 
history acquired in the first scanning position and 
the B-scan, i.e. the waterfall plot collecting the 
time histories recorded at the different angular 
position on the axle circumference, are plotted in 
Fig. 4b and c, respectively. Both the time history 
and the B-scan evidence the Rayleigh wave 
arrival time at about 87 us. However, other waves 
are visible in the waveform which are reflections 
of longitudinal and shear waves due to the 
geometry complexity (section transitions). The B- 
scan close-up around the Rayleigh wave time of 
arrival is shown in Fig. 4d, evidencing the 
presence of the defect, which produces a strong 
attenuation of the Rayleigh wave itself. 

The RMS plot over the time axis (abscissa) of 
the B-scan is shown in Fig. 5. This plot was 
normalized against the maximum RMS value. The 
contrast between the RMS in the damaged and in 
the undamaged areas is 5.1 dB, which was 
calculated considering the minimum values of the 
RMS (in the damaged area) and the floor values 
RMS (in the undamaged area) . 

3. Experimental results on train 
wheel 

The laser ultrasonic technique has been applied 
to a train wheel where Shattered Rim Cracks 
(SRC), e.g. subsurface cracks parallel to the tread 
surface, have been created. Such kind of defect is 
one of the dominant wheel failure types and is the 
result of large fatigue cracks that propagate 
roughly parallel to the wheel tread surface, at a 
depth of about 12-20 mm beneath the thread. In 
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practice, the defect was simulated with a 7 mm 
diameter flat-bottom hole generated perpendicular 
to the rim face, to a depth of 90 mm and with the 
axis at a distance of 10 mm from the inside 
diameter of the rim. 



a. Experimental set-up 



Experimental set-up - defect D4 



AIR-COUPLED rv 
ULTRASOUND |L 




b. Time history close-up on the elastic waves time of arrival 



70 75 80 85 90 95 100 




d. B-scan close-up around the Rayleigh wave time of arrival 

■ 0.04 M 



s 



Fig. 4: Experimental set-up scheme (a), Time 
histories (b), B-scan (c) and (d) 




20 40 60 

Angular position [deg] 

Fig. 5: Normalized RMS on D4 in function of the 
scanning angle 



3.1 Test bench 

The experimental set-up is shown in Fig. 6, for 
both the tested configurations. The pulsed laser 
for the ultrasound generation and the receiving 
probe were the same as the ones exploited for the 
train axle and described in Section 2.1 . 

The laser source was pointed on the outer edge 
of the wheel and the ultrasound probe was located 
on the thread, for the configuration CI , and on the 
inner edge, for the configuration C2, as shown in 
Fig. 6. In both configurations the system was 
working in transmission mode and the presence of 
the defect was detected from its interaction with 
the ultrasound bulk waves. The probe axis was 
tilted with respect to the surface normal of 6 deg. 
That angle was optimized for the detection of the 
transversal wave which resulted as the most 
efficient propagation wave in the wheel test case. 




AIR-COUPLED 
ULTRASOUND 
PROBE 



Configuration C1 




Configuration C2 



Fig. 6: Experimental set-up scheme for two 
measurement configurations 

3.2 Analysis of results 

The results obtained for the two configurations 
analyzed are reported hereafter. Fig. 7 shows the 
time history zoomed in the bulk waves time of 
arrival and the B-scan given in the same time span 
for configuration CI. The longitudinal and 
transversal wave time of arrival are 58 and 70 lis. 
The associated wave trains are visible in both the 
time history and the B-scan. Further propagation 
waves, related to reflection and scattering 
phenomena occurring in a complex geometry 
such as that of the wheel, are visible in the signal. 
The defect acts on the bulk waves as a strong 
attenuator due to the fact that the waves pass 
through a medium with different impedance. 
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The normalized RMS calculated over the time 
axis on a time window of 30 fis, from 60 to 90 (is, 
for all the angular positions is shown in Fig. 7c. 
The contrast between the RMS in the damaged 
and undamaged areas is 8.1 dB. 

The results obtained in the configuration C2 are 
reported in Fig. 8 where the time history, the B- 
scan and the RMS plot are given. The defect is 
well evident in the B-scan and in the RMS plot, 
even though the attenuation is less (5.7 dB) than 
the one produced in the configuration CI (8.1 
dB). However, C2 would be a more suitable 
configuration for practical use. 



Configuration C1 



a. Time history 




| -0.01 



50 60 70 80 90 
Time [us] 



b. B-scan 




< 8 



c. Normalized RMS (30 us, from 60 to 90 us) 




2 4 6 

Angular position [deg] 



Fig. 7: Ultrasound time history, B-scan close-up 
around the transversal wave time of arrival and 
RMS plot in function of the scanning angle - 
Configuration CI 



4. Sensitivity analysis to surface 
conditions 

When NDT methods based on laser ultrasonics 
are applied in field to working structures, the 
surface condition of the test item influences the 
diagnosis and it can produce false positive. 



Configuration C2 



a. Time history 




60 70 80 
Time [us] 




70 80 90 
Time [us] 



c. Normalized RMS (30 us, from 60 to 90 us) 




2 4 6 

Angular position [deg] 



Fig. 8: Ultrasound time history, B-scan close-up 
around the transversal wave time of arrival and 
RMS plot in function of the scanning angle - 
Configuration C2 

This may be particularly critical when Rayleigh 
wave propagation is exploited for damage 
detection, as it was done in the axle investigation 
reported in Section 2. In order to evaluate the 
interaction between the ultrasound propagation 
and the surface condition four different conditions 
have been considered: 

- natural oxidized surface on the laser 
impinging area (thickness of about 35 um); 

- areas with natural oxidation between the 
laser source and the receiving probe (extension 
of about 16 mm in the axial direction); 

- painted surface on the laser impinging area 
(acrylic paint coating with thickness of about 45 

- painted areas between the laser source and 
the receiving probe (acrylic paint coating with 
extension of about 40 mm in the axial 
direction). 

Those different configurations have been tested 
on the train axle described in Section 2 and the 
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experimental set-up is illustrated in Fig. 9a. The 
laser source and the receiving probe were installed 
at a distance of 136 mm. The ultrasound probe 
was fixed at a distance of 34 mm from the axle 
surface. The position of the laser source and the 
receiving probe were kept constant for all the 
configurations. The time histories recorded when 
the surface was clean, when the laser impinged on 
the painted surface and when painted areas were 
existing between the laser source and the probe 
are illustrated in Fig. 9b. It is evident the 
ultrasound is subject to an attenuation when the 
surface condition is not optimal. The attenuation 
of the RMS value is of 3.4 dB when the laser 
impinges on a painted surface and of 0.3 dB when 
the paint is in the surface path between the laser 
source and the probe. 

Fig. 9c shows the time histories recorded when 
the surface was in good condition, when the laser 
was focused on an oxidized area and when the 
oxidation was present in between the laser source 
and the probe. The attenuation is evident again 
and it is of 4.6 dB when the laser interacts with 
the oxidized surface and 1.1 dB when the 
oxidation is in the path between the laser source 
and the probe. 

The analysis brings to conclude the surface 
condition is an important parameter to design an 
optimal laser ultrasonics test, when considering 
the laser impinging area. On the other hand, the 
system is less sensitive to presence of oxidation or 
painted area in the propagation path. 

5. Conclusions 

The paper has shown the applicability of laser- 
ultrasonics, a non-destructive, non-contact 
technique, for the inspection of high speed train 
hollow axles and wheels. The main advantage of 
the proposed technique with respect to the state of 
the art is that it operates in thermo-elastic regime 
and therefore can be considered a non-destructive 
method. Furthermore, it is based on the fact that 
the defect produces an attenuation of the direct 
Rayleigh wave, when working in reflection mode, 
as for the axles, and of the bulk waves, mainly the 
transversal one, when working in transmission 
mode, as for the wheel. That effect, i.e. the 
attenuation of the propagating waves, can be 
therefore used as defect classification feature also 
with very complex geometry. 

A convex defect typically generated by fatigue 
on the external surfaces of an axle has been 
investigated. Such crack has been identified by 
observing the Rayleigh wave propagation, it 
presenting an evident attenuation when passing 



through the damaged area. The contrast between 
the undamaged area and the damaged one reached 
6.5 dB and made it possible to identify the 
presence of the defect. 



a. Experimental set-up 



RECEIVING PROBE 
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Fig. 9: Experimental set-up scheme (a), 
Ultrasound time history for sensitivity analysis to 
the paint (b), Ultrasound time history for 
sensitivity analysis to the oxidation (c) 
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On the wheel a Shattered Rim Crack has been 
created on purpose to reproduce typical 
subsurface defects that propagate roughly parallel 
to the wheel tread surface along the wheel tread. 
This defect acted as strong attenuator of the bulk 
waves and in particular of the transversal one, the 
technique working in thermo-elastic regime. The 
contrast between the damaged and undamaged 
area was between 5.7 and 8.1 dB depending on 
the laser ultrasonics configuration. 

A sensitivity analysis to the surface condition 
has been also performed to consider the influence 
of oxidation and painting on the propagation of 
the Rayleigh wave and on the laser ultrasonics 
results. It has been shown that when the altered 
surface is located in the laser impinging area, the 
induced attenuation is relevant (3 .4 dB for painted 
area, 4.6 dB for oxidized area) and can easily 
produce false positive in the diagnostic procedure. 
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Abstract 

The purpose of the paper is to provide an analysis 
of the attitude exhibited by demodulation 
algorithms of Fiber Bragg Gratings (FBG) based 
on Cross Correlation methods. The Cross 
Correlation technique here considered is based 
on calculations performed in spectral domain. Its 
behavior against the peak-locking phenomenon 
has been analyzed and a way for attenuating it 
has been devised. The new Iterative Cross 
Correlation algorithm improves the FBG 
demodulation accuracy and precision, giving 
remarkably attenuated peak-locking amplitude. 

1. Introduction 

The Fiber Bragg Gratings (FBGs) are optic fiber 
sensors used for local measurement of strain, 
temperature, and/or pressure. Its reflection 
spectrum exhibits a peak at the so-called Bragg 
wavelength, where the fiber shows the maximum 
reflectivity. When a FBG experiences a 
perturbation due to a measurand change, its 
reflection spectrum is subjected to a shift in the 
wavelength domain and a distortion. 

An accurate and precise evaluation of such shift 
allows the calculation of the perturbing 
measurand changes . 

Typically, a FBG spectrum is obtained by means 
of a static optical interrogator which is based on 
spectral reconstruction technique [2] [3] [5] . This 
type of configuration consists of a laser light 
source (characterized by a tunable wavelength), 
swept over the wavelength range (usually 1510 
nm and 1590 nm), with a sampling resolution AX. 
This wavelength sampling resolution limits the 
FBG spectra peak tracking operation [4] [6] . 

The FBG spectrum is achieved once the 
synchronizaztion between the back-reflected light 
and the swept laser source takes place. 

Once the perturbed FBG spectra are acquired, a 
method for tracking the spectrum wavelength 
shift (which allows to determine the measurand 
change) is required. Typically such task is 
performed in the wavelength domain, and the 
spectrum wavelength shift is detected by means of 
a least-squares fit of the discretized spectrum in 



such a way that an accurate and precise estimation 
of the Bragg wavelength shift occurring between 
two FBG spectra is achievable [6] [7] [10]. Such 
technique, among the others operating in the 
wavelength domain, are quite computationally 
efficient and lead to acceptable accuracies and 
precisions. 

2. Cross Correlation Method (CC) for 
FBG demodulation based on 
spectral techniques 

In this paper, a novel technique for FBG 
spectrum shift tracking based on Cross 
Correlation technique is devised. 

There are several papers in the scientific 
literature dealing with and developing such 
tracking techniques [1] [8] [9]. 

The time consumption performance granted by 
the Cross Correlation method is comparable to the 
one of the standard methods operating in the 
wavelength domain and its strength is given by 
the remarkable performance in terms of accuracy 
and precision in the FBG wavelength shift 
detection [8] . 

In order to reduce the computational complexity 
characterizing such technique, the Cross 
Correlation calculation considered takes 
advantage of the Cross Correlation Theorem and 
the Convolution Theorem. They state that the 
Fourier Transform of the cross-correlation 
between two signals defined in their own physical 
domain equals the product of the individual 
Fourier Transforms of both signals, one of which 
is complex-conjugated [8] . Given N the length of 
the signals, the overall computational complexity 
is OfNlogN] instead of 0[N 2 ], typical of Cross - 
Correlation calculations based on the standard 
moving sliding-dot product. 

The mathematical model used to represent the 
input FBG spectra (the reference spectrum and the 
shifted one because of an applied perturbation) is 
given by the following relations 

x(A) = pW + r^iA) 

Eq. 1 

y(X) = p(A + AX) + n 2 (A) 
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In the previous relations, the term AX refers to 
the Bragg wavelength shift, which allows the 
calculation of the perturbing measurand. The 
terms x(X) and y(X) refer to the reference and 
shifted FBG spectra, respectively. The term p is 
the un-noised component of the FBG spectra 
(reference and shifted ones), while the terms ni(X) 
and m(X) refer to the noise components of both 
spectra. 

The Cross Correlation function developed in [8] , 
and computed in the Fourier domain, is given by 
Eq.2. 

R xy (A) = T- l i [ T[y(A)].T[x(A)J ] j Eq.2 

In Eq.2 the Fourier and the Inverse Fourier 
Transform operators have been introduced. 




1539.8 1540 1540.2 1540.4 1540.6 1540.8 1541 

Wavelength [nm] 



Fig. 1: Example of FBG reference and shifted 
spectra for a Bragg wavelength shift of ~114 pm. 




0 50 100 150 200 250 300 350 



Bragg Wavelength Shift A^ B [pm] 

Fig. 2: Cross Correlation function of FBG spectra 
in fig. 1. The pictured points are subjected to the 
Gaussian fit for the determination of interpolated 
Cross Correlation peak x-position. 



Correlation function peak. Due to Cross 
Correlation function discrete nature, the peak 
interpolation carried out by means of 3 -point 
Gaussian fit has been performed. The authors 
propose the implementation of such operation, in 
order to get better accuracies and avoid the FBG 
spectral resolution limitations. 

2.1. Cross Correlation Performance 
Simulation through Monte Carlo 
Method 

In order to evaluate the enhancements led by the 
described technique (mainly in terms of better 
accuracies and precisions) the authors, in a 
previous paper [8] , applied the procedure to FBG 
spectra simulated by means of the Transfer 
Matrix Method, a powerful mathematical tool 
allowing the simulation of the FBG behavior 
under several strain conditions. 

The physical properties of the FBG here 
considered are summarized in the following table. 
Table I 

Physical parameters of the simulated uniform fbg 



Symbol Description Quantity 

L FBG total length 10 mm 

ricore core refractive index 1 .46 

n c iadding cladding refractive index 1 .44 

r core core radius 4 um 

r •cladding cladding radius 62.5 um 

n e ff effective refractive index 1 .46 

u \ i *• • * 0.26,0.11: 

p e ,pii,pi2 pnotoelastic coefficients 0 252 

v Poisson ratio 0.17 

AX spectral resolution 10 pm 

so applied strain coefficient 10~ 4 us 



The simulation produced a set of spectra which 
exhibit Bragg wavelength shifts from 2.34 pm up 
to 114 pm, with a FBG spectra resolution of 10 
pm. 

Fig. 3 shows the reflection spectra as the FBG is 
strained. The strain function considered is 
linearly -time varying. 



Fig. 2 shows the operations performed on the 
Cross Correlation function, in order to achieve the 
Bragg wavelength shift that is given by the Cross 
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Reflection Spectra of FBG under Strain 




1540.6 



1540.7 



1540.1 1540.2 1540.3 1540.4 1540.5 

Wavelength [pm] 

Fig .3 Some simulated reflection spectra for the fiber 
under consideration. The fiber considered 
experiences a linearly time varying strain. 

In order to evaluate the Cross Correlation 
performance in noisy conditions (since the 
acquired spectra are actually affected by noise), 
the Monte Carlo method have been implemented. 
It allows the simulation of the FBG strain-shifted 
reflection spectra for several noise levels, 
expressed in terms of spectrum Signal-to-Noise 
Ratio (SNR). Thus, each shifted spectrum is 
considered white Gaussian noised with a SNR 
from 20 dB up to 60 dB with a 2 dB step. The 
Monte Carlo implementation has the purpose of 
evaluating each algorithm performances for 
several spectral noise levels. 



Number of Repetitions: 1000 



FBG Spectrum 
Generation by 
means of TMM 



White Random 
Gaussian Noise 
Generation by 
means of AWGN 
Algorithm 



Implementation of 
FBG demodulation 
► Algorithm and 
Estimation of Bragg 
Wavelength Shift AX 



Fig. 4 How the implemented Monte Carlo method 
works. 

The paper [8] introduces a series of comparisons 
between the described Cross Correlation 
technique (see Eq. 2) and the usually employed 
FBG demodulation techniques operating in the 
wavelength domain (in particular, the paper 
considers the quadratic least-squares fit and the 
centroid algorithm, introduced by [4]). 

The carried out performance analyses involve 
the time consumption, the accuracy and the 
precision indices. The accuracy and the precision 
indices are given by the following relations. 



i M 

RMSE = A — V ( A>i 



calculated effective 



) 

Eq.3 



1 M / 



Eq.4 



The Root Mean Squared Error and the Standard 
Deviation values are calculated for each sample 
(containing M=1000 simulated spectra) and for 
each SNR considered. 

Figg.5 and 6 show the RMSE comparison 
between the Cross Correlation method described 
and the other standard FBG demodulation 
techniques. 

Root Mean Squared Difference for input Spectra SNR 20 dB 



- Quadratic Least Squares Fitting 

- Centroid 

- Cross Correlation Method 




Bragg Wavelength ShiftAA. B [pm] 

Fig. 5 RMSE trend for the three methods vs the 
Effective Bragg wavelength shift. 



Root Mean Squared Difference for AX B =114 pm vs SNR 




Quadratic Least Squs 
Centroid 
Cross Correlation Method 



Input Spectra SNR [dB] 

Fig. 6 RMSE comparisons for an effective Bragg 
wavelength shift ~114 pm vs spectra SNRs. 

The pictures show that for a given SNR the 
Cross Correlation algorithm perform better and 
for a given Bragg wavelength shift the RMSE 
granted by the Cross Correlation method is 
always less than the other two methods, 
regardless of the noise. 

The following table shows a computational time 
comparison. These comparisons have been 
performed on an 8 Gb memory Intel® Core i7™ 
3632QM (2.20 GHz) computer. 
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TABLE II 

Av erage computation time required by each met hod 
Method Average Time Required 



Quadratic 
Least- 
Squares Fit 



1.5-10" 3 s 



Centroid 
Algorithm 



7.2-10" 4 s 



Proposed 

Spectral 

Cross 

Correlation 



1.6- HTs 



As can be noted, the Cross Correlation method 
offers also the best time consumption 
performance. 



s Correlation FBG spectra 20 dB 




10 20 30 40 50 60 70 



100 110 120 



Effective Bragg Wavelength ShiftAX- B [pm] 

Fig. 7 RMSE trend for Cross Correlation method vs 
effective Bragg wavelength shift for SNR=20 dB 
revealing the peak-locking phenomenon. 

Despite the excellent performance typical of the 
Cross Correlation method, the oscillatory trend 
exhibited by the RMSE index (as the effective 
Bragg wavelength changes) could lead to 
inaccurate Ahmgg evaluation. The minimum error 
is located at integer multiples of sampling spectral 
resolution (10 pm in this case). The calculated 
AlBmgg then drifts away from the actual one 
following an oscillatory trend with a period 
equaling the spectral resolution. 

This particular trend resembles the one occurring 
in Particle Image Velocimetry (PIV), as described 
by [11] [12]. 

This phenomenon is due to the inability showed 
by the FBG demodulation technique, to resolve 
the effective Ahmgg whenever it is less than one 
sampling spectral resolution. This inability is 
accentuated by larger sampling spectral 
resolutions and is influenced by the 3 -point peak 
fit of the Cross Correlation function (with a 
Gaussian fit as considered in this paper). Thus, 
the peak-locking effect (definition actually 
borrowed from PIV) can be considered as a bias 



towards integer values of sampling spectral 
resolution. 

3. Iterative Cross Correlation Method 
ICC) 

In order to reduce the AXsmgg evaluation 
inaccuracies due to the peak-locking effect, the 
authors investigate the enhancement of the 
described Cross Correlation algorithm, by 
performing it iteratively. 

The proposed improvement is going to be 
discussed later. 

3.1. Further Considerations on Peak- 
Locking Effect 

This paragraph is aimed to provide further 
insight about peak-locking effect. 

As previously said, the accuracy parameter 
RMSE reveals a quasi-periodic trend that should 
be considered as a biasing phenomenon induced 
by the sampling spectral resolution and the 
interpolation operations based on peak-fit 
performed on the Cross Correlation function. 

Such phenomenon is due to the relative inability 
of the demodulation techniques based on FFT- 
based Cross Correlation function, to correctly 
detect the Bragg wavelength shifts at sub- 
resolution extent [8] 

Indeed, the CC function has the resolution of the 
input spectra. Therefore, the ability of accurately 
finding the CC function peak is limited by its own 
resolution. This inability is partially overcame by 
the 3 -point peak fit operation. 

Furthermore, if the defined effective AXsmgg is an 
entire multiple of resolution, the 3-point Cross 
Correlation peak fit returns a peak point whose 
abscissa remains mostly constant (as the added 
random Gaussian noise superimposed to both 
spectra changes during the Monte Carlo 
iterations) and it accurately represents the 
effective Bragg wavelength. 



c 12.8 
I 

£ 12.6 

| 12.4 

o 

1 122 
1: 



Fig. 8 Cross Correlation peak interpolation for an 
effective AXBmgg = 100 pm. Here the interpolated 
peak is coincident to the peak value of discrete 
Cross Correlation function 
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Fig. 9 shows the peak fit process when the 
effective AXsmgg is not equal to an entire multiple 
of spectral resolution (105 pm in this case) 

14i 1 1 1 1 1 1 1 1 1 

| 13.5- 

1 12.5 - 

^ I 1 

o 115 _ Gaussian Fit 

i 

O 11 - 

1 °'80 85 90 95 100 105 110 115 120 125 

Effective A*. Bragg [pm] 

Fig. 9 Cross Correlation peak interpolation for an 
effective JlBmgg = 105 pm. The star indicates the 
peak of the Gaussian fit. 

As the Monte Carlo iterations evolve, the peak 
of the discretized Cross Correlation function and 
the peak of the gaussian fit are not coincident and 
their mutual position randomly changes during 
the Monte Carlo iterations. This behavior is 
mainly responsible for the RMSE relatively high 
values for effective AXsragg not equal to a multiple 
integer of resolution. 

In order to provide a quantitative estimator of 
such phenomenon, the authors implement a series 
of statistical considerations on the already defined 
RMSE accuracy index. 

The overall accuracy of calculated Bragg 
wavelength displacement is a combination of 
various aspects, and it can be explained (for a 
specified SNR) by means of the following 
relation. 

A/l. = AX eff + P + £. Eq.5 

In the Eq. 5, the term Ah is each of the 
calculated Bragg wavelength shift in each sample 
generated by Monte Carlo simulation for a 
constant SNR (sample that contains M=1000 
simulated FBG spectra), AX e jj is the effective 
Bragg wavelength shift, the term ft refers to the 
bias error and Si is the random error with zero 
mean. Each sample (for each effective Bragg 
wavelength displacement and SNR and containing 
M generated spectra) is such that the average of 
the random error is null and the bias error, /?, is 
constant. 

Therefore, the average value of the calculated 
Bragg wavelength shift for each Monte Carlo 
generated sample can be calculated as follows. 

Y m ^ M 

AX = — VAA =AA~ + fi + — Ye. = AA ff + 8 

Eq. 6 



The random error level is evaluable by means of 
the root mean squared fluctuation of the 
calculated values around the mean value. 

Pi ^ 2 
Let us recall the RMSE value (Eq.3). 

r~j m 2 r~j m 

*MSE = ^-g(M-A^) =^MW + ei) 

By means of Eq. 8, the bias term ft, 
representative of the peak-locking phenomenon 
can be calculated. 

3.2. Description of the ICC Method 

In order to attenuate the oscillatory trend 
exhibited by the RMSE index as the AXsmgg varies, 
the author propose an iteratively performed Cross 
Correlation method. 

The peak position (whose abscissa indicates the 
computed Ahragg) is up to date at each iteration. 
The peak position detected in the latter iteration is 
used as one of the three-interpolation point taken 
close to the peak position detected in the current 
loop. 

At the first step, a 3 -point peak gaussian fit of 
the discretized Cross Correlation function is 
performed. Once the peak of the interpolating 
function is calculated, its x- value is used for the 
calculation of the renewed Cross Correlation peak 
amplitude from the Fourier coefficients. 

The iterations stop until the error function 
(calculated as the absolute value of the difference 
between the last two calculated Cross Correlation 
peak amplitudes) reaches a desired value. 

This method is more time consuming depending 
on the desired accuracy level, but it allows the 
achievement of better accuracy and precision 
performance. 

The following flow chart shows the several 
calculations involved in the proposed method. 
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( START ) 



Cross-Correlation calculation and its peak 
determination 



3 -point Gaussian peak fit 

1 



Recalculation of peak position and height and 
adaptive adjustment of peak boundaries 



Calculation of error function as difference 
between last peak position and current one 




Fig. 10 Proposed Iterative Cross Correlation 
method flow chart. 



4. Performance Comparisons 
between proposed ICC and the CC 
Methods 

The following picture shows the bias error trend 
against the effective Bragg wavelength shift for 
CC and ICC methods. 



— Cross Correlation FBG spectra 30 dl 
Iterative Cross Correlation 30 dB 

— Cross Correlation 60 dB 

— Iterative Cross Correlation 60 dB 




10 20 30 40 50 



100 110 120 



Effective Bragg Wavelength ShiftA^ B [pm] 

Fig. 11: Bias comparison between the Cross 
Correlation and the Iterative Cross Correlation 
methods for input FBG spectra with SNR 30 dB 
and 60 dB. 

By means of an iteratively performed Cross 
Correlation algorithm, a significant reduction of 
the oscillatory trend is obtained. This 



improvement is more evident for low and very 
low SNRs, where the best oscillatory trend 
reduction could reach values up to 80%. 

The following table reports the average peak-to- 
peak values of the bias error index (whose trend 
is considered a waveform) for both methods under 
comparison and for several SNRs. 

Table III 



Method 


20 dB 


30 dB 


60 dB 


CC 


0.18 


0.04 


0.07 


ICC 


0.11 


0.01 


0.006 



The table, as Fig. 11, shows a significant 
decrease of the peak-to-peak value of the 
waveform associated to the bias error trend as the 
noise level decreases. The efficacy of the 
proposed Iterative Cross Correlation method 
become, thus, more robust with high FBG spectra 
SNRs. 

5. Some Comments about Peak- 
Locking in PIV 

In the PIV technique, several physical 
parameters contribute to the peak-locking effect. 
In this case, the average seeding particle diameter 
plays a crucial role. Indeed, when the seeding 
particle diameter becomes significantly small, the 
detected displacement between two time- 
consecutive frames is biased towards integer 
values (especially when the particles are smaller 
than the interrogation pixel dimension). 

The effect exacerbates when a gaussian sub- 
pixel peak estimator is employed. 

In the case of 1-D signals, (such as the ones 
contemplated in this paper and relative to the 
FBG reflection spectra) the 3 -point gaussian peak 
fit has an improving effect (since it allows 
overcoming the limitation due to not fine 
sampling spectral resolutions). In 2-D signal, in 
imagery and PIV imagery, the concomitance of 
pixel discretization, limited by seeding particle 
diameters, seeding particles density and gaussian 
peak fit accentuate the peak-locking trend. 

In addition, the biasing phenomenon 
characterizing the PIV techniques is exacerbated 
by the FFT-based Cross Correlation calculation, 
which can be described by the following block 
diagram. 
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1 st Input ID (FBG 

spectra) 
1 st Input 2D (PIV 

images) 



2nd Input ID 
(FBG spectra) 
2nd Input 2D (PIV 

images) 



Real-to-Complex 
FFT 



Real-to-Complex 
FFT 



Complex Cojugate 
Multiplication 



Complex-to-Real 
Inverse FFT 



Cross Correlation 
Data 



Fig. 12: Block Diagram relative to CC calculation 
for both ID signals (i.e. FBG spectra) and 2D PIV 
images. 

Although the CC computation is relatively easy 
to perform in the Fourier domain (with the main 
advantage of significantly reducing computational 
complexity), the pixel discretization in PIV filed 
and the resolution limitation for FBG signals, in 
addition to finite interrogation window length, 
introduces errors within the estimation of the 
spectra content, inhibiting the performance of the 
Fourier-based calculations. Furthermore, the 
periodicity of the assumed signal boundaries 
implies aliasing effects, affecting negatively the 
ability of the demodulation method to accurately 
detect shifts not equal to integer multiples . 

In PIV applications two possible stages of 
operation have been introduced, [11] [12]. The 
first one deals with the pre-processing of the 
images under analysis. The operations involved in 
this step are image windowing and zero-padding. 
The higher level of operation involves an 
iteratively performed CC calculation on input 
images that are manipulated in order to obtain (at 
the end of each iteration) successive 
displacements tending to sub-pixel levels. Thus, 
in the PIV technique, the joint implementation of 
an iteratively performed CC calculation and an 
adequate manipulation of the input images under 
analysis, successfully leads to a sensible 
attenuation of the bias error characterizing the 
peak-locking effect. 



6. Conclusion 

In this paper, the authors set the purpose of 
enhancing the Cross Correlation function 
(addressed to FBG signals demodulation for 
determining AXBragg), which is calculated by 
means of the Fast Fourier Transform tool, 
according to the Convolution Theorem. In 
particular, the proposed enhancement has the aim 
of reducing the so-called peak-locking effect. This 
term, borrowed from the Particle Image 
Velocimetry, refers to a biasing phenomenon, 
which leads to Bragg wavelength shift evaluations 
inaccuracies. 

The physical and mathematical causes implying 
such behavior have been explained, for the FBG 
spectra. 

The authors show that an iterative execution of 
the Cross Correlation algorithm (in conjunction 
with the determination of AX B rag g at sub-resolution 
extent by means of 3 -point peak gaussian fit) 
leads to more accurate AXBragg evaluations, 
especially for effective Bragg wavelength shifts 
not equal to entire multiples of the FBG sampling 
spectral resolution. 

Several FBG reflection spectra have been 
generated by means of the Transfer Matrix 
Method, in order to simulate the fiber behavior 
under an applied strain. Therefore, in order to take 
into account the noise (inevitably characterizing 
the effective sampled spectra), some white 
gaussian noise have been added and the efficacy 
of the proposed algorithm is then tested on a 
series of FBG spectra, generated by means of 
Monte Carlo method. 

The overall results, in terms of accuracy and 
precision of the detected AX Br agg, have been 
compared with those relative to the Cross 
Correlation technique applied once. 

The comparisons between the ICC and the CC 
show a remarkable improvement of the accuracy, 
precision and a significant reduction of the 
fluctuations exhibited by the accuracy index, 
symptom of the peak-locking effect. 
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Abstract 

Focus of the present work is the evaluation of the 
air flow rate in a pipe of squared cross section 
through the direct measurement of the mean 
velocity at a fixed location of the section. To this 
aim, an extension of the methodology already 
utilized for circular section pipes is here 
proposed. The method simplifies and improves the 
procedure suggested by the technical standard EN 
12599. Specifically, it has been possible to state a 
geometrical position in the cross section where to 
straight measure the mean flow velocity in order 
to evaluate the flow rate. Here, it is shown that 
the proposed methodology allows the estimation 
of the flow rate with higher accuracy than the 
normed one does, estimating errors less than 
those accepted in common practice. 

1. Introduction 

The measure of the flow rate is very important in 
scientific, industrial and economic fields; in 
particular, the measure of the air flow rate and the 
evaluation of the leakages in piping networks 
could be important for the energy efficiency of 
the complete air conditioning system when trying 
to improve energy conservation. For several 
constructional reasons, this type of installations 
have, normally, a rectangular/squared sections. 

The most straightforward method for measuring 
the flow rate is based on its definition and consists 
in measuring the volume or mass of the fluid 
flowing in steady conditions in a given time 
interval. This method is very time consuming and 
laborious (it has in fact no industrial applications; 
is used, for example, in the process of primary 
calibration of instruments), so the flow rate must 
be evaluated by measuring the local velocity (or 
the average one), knowing the appropriate cross- 
section. Usually, the measure of local velocity of 
the fluid in pipes is obtained by means of Pitot 
Static Tube, wire anemometer, and so on, while 
the measure of the average velocity of the flow is 
obtained by the using throttling devices, 
rotameter, electromagnetic flow meter, turbine 
meter, ultrasonic meter, etc. 



In literature there are many papers dealing with 
the problem of measuring the flow of fluids [1]- 
[4]; in [5] it is experimentally shown that, for 
turbulent flows and circular ducts, the ratio 
between the punctual speed measured at 
approximately 3 A of the radius of the duct from its 
axis and the average speed in the section is, at 
different values of the Reynolds number, 
approximately unitary. This leads to the 
conclusion that it would be enough to measure the 
velocity right at a point in an opportune cross - 
section of the pipe (punctual speed) to obtain 
from that the average speed, and then the flow 
rate. Starting from these considerations, the 
authors derive empirically a similar relationship 
between the punctual speed and the average speed 
of a turbulent flow in ducts with squared cross- 
section exploiting the concept of equivalent 
circular cross section. 

In this work experimental results are presented 
relatively to flow rate measurements in a square 
cross-section tube; in particular, the authors 
evaluate the flow rates from velocity 
measurements (acquired by LDA instrumentation) 
according to the EN 12599 Standard guidelines 
and derive an empirical formula suitable to obtain 
the position in the section where directly 
measuring the flow rate velocity. 

The results obtained show that the proposed 
method provides a low measurement uncertainty 
for flow rate evaluations in turbulent regime and 
may be applied to air conditioning systems as 
well as, more generally, to installations using 
ducts of squared cross-section. 

2. Empirical model 

2.1. Background: Circular cross 
section ducts 

It is know that for all continuous fluids the 
velocity at the wall of the duct is stationary 
relative to it (no slip condition). 

The velocity then increases by moving toward 
the axis of the duct, both for laminar and turbulent 
flow, Fig. 1. 
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In laminar case and for fully developed profile, 
the mean velocity of the flow in a circular cross 
section pipe is given by 

Vo 
2 

where R is the cross section radius and 

lApR 2 
Vo=- 



4 



Eq.2 



If the flow is turbulent, Eq.l is no more valid 
and it has been proved that the local flow velocity 
can be evaluated by the empirical formula 



Vn V RJ 



Eq. 3 



where n is linked to the Reynolds number from 
experimental data [5] and here reported in Table 1 

Table 1 Reynolds number versus n 



Re 4.0 10 3 


2.3 10 4 


1,1 10 5 


1.1 10 6 2.0 10 6 to3.2 10 6 


n 6.0 


6.6 


7.0 


8.8 10.0 



Fig. 1: Laminar and turbulent pipe velocity profiles 

In such condition, the mean flow velocity is 
measurable at about 3/4 of R from the axis (i.e., 
see Preston device) and then, 
3 V 

r = -R -» 1 EqA 
A V H 

Integrating Eq.3, the mean velocity can be 
expressed by: 
V _ 2n 2 

Vb"(n + l)(2n + l) Eq - S 
or 

l 



V (n + l)(2n + l) 



6 



7 2n 2 

Combining these last two equations, the velocity 
profile descends as given by Eq. 3 . 

Experimental data show that V = V at r = 
0.758/? for turbulent flow. 



2.2. Squared cross-section ducts 

In the case of turbulent flows, it is possible to 
extend the previous model to squared cross 
section pipes in the form: 



l 

V ( 2x\n 



Eq.l 



being x the distance from the axis of a generic 
point in the pipe cross-section and I is the side of 
the squared section considered. This empirical 
formula in some way descends from the Eq. 3, by 
suitably rearranging the quantities there present 
for similarity. 

By considering the equivalent circular section 
draining the same volumetric flow rate, such as: 
VS r = VS eqc Eq.8 
where S r and S eq c are the squared and the 
equivalent circular cross-sections respectively, the 
radius of the latter is given by 



Rpn r 



Eq.9 



where / is the side of the squared cross-section. 
The direct proportionality between R eq _ c and / 
clarifies the proposed formula where the 
coefficient 2 approximates the square root of n. 
Eq. 9 then allows to apply the original model for 
the equivalent circular cross -section to the actual 
squared one and to find the appropriate location x 
in the squared section area where to measure the 
expected mean velocity for the flow rate 
evaluation. The use of Vn instead of 2 multiplier 
in Eq. 7 gives more accurate results, of course, 
but the discrepancy is very insignificant in the 
frame of the practically used uncertainties. 

3. Experimental tests 

The experimental tests have been carry out to 
prove the reliability of the proposed model for 
locating the right position in a squared section 
where to measure the mean velocity and so 
evaluating the volumetric flow rate flowing inside 
the duct. 

3.1. Experimental setup and 
instrumentation 

A small subsonic wind tunnel has been arranged 
whose test section has been properly designed, 
Fig. 2. 

The tunnel is an open circuit, in which a 
maximum air velocity of 32 m/s has been 
measured along the axis of the test section. 
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Test section 











40 I 



Fig. 2: Test section in the subsonic wind tunnel 

The inlet nozzle, 1.125 m long, has a contraction 
ratio of 4:1 (in terms of section sides ratio), while 
the test volume is 0.600 m long having a squared 
transversal section with side of 0.153 m. Its walls 
have been realized by plexiglass except that one 
crossed by the laser beam, a glass plate 4 mm 
thick, in order to limit diffraction phenomena 
across it. The diffuser, 1.200 m long, has a 
squared inlet section and a circular outlet section 
(<|> = 0.355 m). At the exit of the diffuser an axial 
fan is connected driven by a three-phase 
asynchronous electrical engine (2900 rpm 
nominal speed), whose round speed is controlled 
by inverter. The test section has been chosen at 
0.400 m from the nozzle exit with the optical 
access for the LDA device used to measure the 
flow velocity profile on it, Fig. 3. 




Fig. 3: Wind tunnel optical access for the LDA 
device. 



3.1.1. LDA and moving system 

A LDA anemometer, FlowLite 2D System by 
Dantec Dynamics, [6] (Fig. 4) has been employed 
having blue (Laser Sapphire, X = 488 nm) and 
green (Laser Yag, X = 532 nm) beams 
combinations, for a power of 200 mW and 
2.2 mm beam diameter. 




Fig. 4: Flowlite LDA System 

The laser probe moves along two binaries 
allowing 4.5 + 0.1 mm horizontal and/or vertical 
displacement steps. The probe optical lens has 
<p = 60 mm and / = 160 mm. The signal 
processor of the LDA is the BSA F60 - Dantec 
Dynamics. 

3.1.2. Seeding generators 

The pneumatic seeding generator here utilized 
was fed by air extracted downstream of a pressure 
controller and located upstream of the inlet 
section of the nozzle, in axial position. 

The seeding was constituted by a 3% solution of 
ethylene and paraffinic oil in order to guarantee a 
better diffraction of light from the seeding 
particles and so a better SNR. Its flow rate was 
regulated through the air pressure feeding the 
generator. 

3.2. The EN 12599 procedure 

A short description of the normed procedure 
regarding the evaluation of the flow rate of a 
fluid, is described as reported in EN 12599. The 
efficiency of the proposed procedure is attested by 
the computation of discrepancies between the 
measured velocities and their theoretical values 
coming from the model. 

In the technical standard EN 12599:2012 [7], 
titled "Ventilation for buildings - Test procedures 
and measurement methods to hand over air 
conditioning and ventilation systems", is 
illustrated the method for evaluating the air flow 
rate in rectangular ducts (see Annex D), and, 
specifically, the measuring methods and devices. 

In Annex D is specified that the measurement 
cross-section of rectangular ducts should be 
divided into i elements of Sj equivalent area as 
shown in Fig. 5 
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Fig. 5: Areas of equal size where pi is the measuring 
point of the area S t 

The mean velocity is obtained as the arithmetic 
mean of the punctual velocity measured at the 
centroids of the sub-areas Si; but there is also 
written that "the number of measurement point 
depends not only on the geometrical size of the 
cross section but also and decisively on the 
velocity profile. In case of large velocity 
differences, the distance between measurement 
points should be smaller and the measurements 
should be appropriately differently evaluated'. 

According to this advice, the authors have 
divided the squared cross-section under 
investigation in analogy to the partitions of 
circular sections and fixed the points of measure, 
as illustrated in Fig. 6, as example. 

m 




Fig. 6: partition of the test section in areas of equal 
size with indication of the measuring points. 




33. Test procedure 

The measuring points are numbered as shown in 
Fig. 8 so indicating the path followed by the 
moving system of the laser probe, and whose 
coordinates are reported in Table 2. 



y /N 



i 

0- 



O — >- 



16 OO 



15 



1 

-Q- 



Q Q C) 



— ee 



14 13 12 11 

O O Q 



10 



17 18 19 

q22< 21 q 



23 



-e — > — e 



24 



25 



Fig. 8: path of the probe into the test section 

Table 2 Cartesian coordinates related to the 
measurement points, pi 



| Pi 1 2 3 4 5 


| -67.5 


0 


63 


54 


0 


| 67.5 


67.5 


67.5 


54 


54 


6 


7 


8 


9 


10 


| -54 


-22.5 


0 


22.5 


67.5 


54 


22.5 


22.5 


22.5 


0 


l p< l ii 


12 


13 


14 


15 


54 


22.5 


0 


-22.5 


-54 


0 


0 


0 


0 


0 


mmm 


17 


18 


MEM 


mm 


| -67.5 


-22.5 


0 


22.5 


54 


0 


-22.5 


-22.5 


-22.5 


-54 


Pi 21 


22 


23 




25 


0 


-54 


-67.5 


0 


63 


| -54 


-54 


-67.5 


-67.5 


-67.5 



The measurements have been carried out for 
three different values of flow rate (3 different 
revolution frequencies controlled by the inverter, 
10, 25 and 40 Hz); further, at each point five 
measures of velocity have been taken every 60 
minutes and each of these measurements has been 
replicated for 5 consecutive days. For each 
acquisition the number of samples has been about 
N=1000 and the acquisition time T s equal to 10 s. 
The acquisitions have been carried out in burst 
mode. 



Fig. 7: Perspective view of the control volume 
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4. Results and comments 

4.1. Velocity profiles 

The averaged velocity profiles are reported 
below for the 3 different revolution regimes, 
obtained by taking the mean values of those 
measured at the correspondent points on the four 
reference lines (2 axes and 2 diagonals) of the 
cross section. 



5,50 



5,10 



4,90 
4,70 



4,50 



velocity profile 10 Hz 



-100,00 -50,00 0,00 50,00 

Test point location (mm) 



100,00 



Fig. 9: avareged velocity profile obtained at 10 Hz 
revolution regime 



16,00 



15,60 



14,80 



14,00 



Velocity profile 25 Hz 



-100,00 -50,00 0,00 50,00 

Test point location (mm) 



100,00 



Fig. 10: avareged velocity profile obtained at 25 Hz 
revolution regime 



26,00 
25,50 



25,00 
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23,00 



Velocity profile 40 Hz 



-100,00 



-50,00 0,00 50,00 

Test point location (mm) 



100,00 



Fig. 11: avareged velocity profile obtained at 40 Hz 
revolution regime 



4.2. Theoretical velocity profiles 
evaluation 

The test section is squared with side I = 
153,0 mm, which corresponds to a radius of the 
equivalent circular cross-section, by using Eq. 9, 
equal to: 



Rpn r 



I 2 



= 86.3 mm 



Calculating the theoretical profile for the case 10 
Hz by Eq. 7 with values of n as reported in the 
text, [5], and here in Tab. 1, it appears as in Fig. 
12 where it is compared with the experimental 
one. 

This result leads to the consideration that n 
might depend not only on Raynolds number but 
also on the geometric form of the pipe cross- 
section. 

theoretical profile vs experimental profile, 10 Hz 





o,uu 
























5,50 














9r± 


























S 4,50 
















4,00 







-100,00 -50,00 0,00 50,00 100,00 

Test point location (mm) 
• experimental profile • theoretical profile 

Fig. 12: experimental averaged velocity profile 
compared with theoretical profile obtained using 
the values of n as reported in [5]. 

Then, new values of n have been calculated by 

Eq. 3 rewritten as: 

V(r)\ If x \ 
-Tr L ) = -ln[l- — ) Eq. 10 

and using the mean experimental velocity profile 
values for V(r), for each revolution regime: the 
slope of the straight line, then, obtained from a 
least squares regression is the best estimation of 
1 / n, and therefore for n. 
The values of n obtained are reported in Table 3 



In 



Table 3 new value of n 



56943 169535 278957 
457 1317 2740 

where the Reynolds number is calculated from 
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pV r , 



•an^eq 



Eq.ll 

and p = 1.206^, u = 18.24 ■ 1(T 6 P a ■ s, for 
T env = 20,0°C and p env = 100 kPa. 

The mean theoretical velocity profiles for each 
flow rate investigated have been then estimated 
by Eq. 7 by means of the new values of n and 
plotted in Fig. 13, Fig. 14 and Fig. 15, where they 
are compared with the experimental ones. 
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Fig. 13: comparison between experimental and 
theoretical averaged velocity profile at 10 Hz. 
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Fig. 14: comparison between experimental and 
theoretical averaged velocity profile at 25 Hz. 
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As shown in the graphs, the theoretical and 
experimental profiles match very well. 

For each one of the revolution regimes, the 
experimental mean velocity is obtained as the 
average value of the 25 test points ones, as: 

25 

1 v 



25 



2/< 



Eq. 12 



where Vj's are the mean weighted experimental 
values of the local velocities obtained by the 
measurements. Its standard deviation is simply 
calculated by the law of propagation of 
uncertainty, SV mean . 

The mean velocity values with their standard 
deviations for each regime are reported in Table 4 

Table 4 flow rate velocities 





Vjnean.exp 


4.95 


14.73 


24.24 


SV mean (m/5) 


0.07 


0.11 


0.18 



Thus, the location x correspondent to the mean 
velocity of the flux can be now estimated from 
Eq.7. 

For all the 3 fluxes experimented the ratio 
has resulted equal to -0.776, compared 
with -0.758 for circular cross-section. 

4.3. Evaluation of measurement 
uncertainty of the proposed method 

The followings Table 5 contains the deviations 
between experimental velocities and theoretical 
ones for each flow rate. 

The maximum relative deviation estimated is 
0.01%. 

4.3.1. Theoretical mean velocity evaluation 
at 0.776 1/2 

The theoretical mean velocity has been 
estimated, for each one of the rotational regimes 
considered, from the empirical model stated 
(Eq. 7) and compared with the experimental 
correspondent value calculated by means of 
Eq. 12, as shown in Table 6. There, the 
discrepancies between the homologous velocity 
values and relative errors are also present. 

The error estimated appears to be less than 3%, 
so validating the hypothesis of obtaining the flow 
velocity for the estimation of the flow rate 
directly from a unique measurement at a chosen 
position even in the squared cross-section pipe as 
well as in circular cross-section ones. 



Fig. 15: comparison between experimental and 
theoretical averaged velocity profile at 40 Hz. 
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Table 5 Comparison between experimental and theoretical velocity averaged profiles 

10 Hz 





V <. , rm/<?i 

v experimental L"V J J 


v theoretical Y m / S \ 


AT/ \m /<:1 


AT/ IV . , % 
L * v / v experimental /u 


-67.50 


5.074 


5.066 


0.008 


0.2 


-54.00 


5.096 


5.076 


0.020 


0.4 


-22.50 


5.095 


5.086 


0.009 


0.2 


0.00 


5.090 


5.090 


0.000 


0.0 


22.50 


5.090 


5.086 


0.004 


0.1 


54.00 


5.075 


5.088 


-0.008 


-0.2 


67.50 


5.076 


5.066 


0.010 


0.2 


25 Hz 


Test point 
locations [mm] 


^experimental [m/s] 


y theoretical [m/s] 


AT/ [m/s] 


/^experimental % 


-67.50 


15.023 


15.044 


-0.021 


-0.1 


-54.00 


15.083 


15.054 


0.029 


0.2 


-22.50 


15.075 


15.064 


0.009 


0.1 


0.00 


15.060 


15.068 


-0.008 


-0.1 


22.50 


15.072 


15.057 


0.015 


0.1 


54.00 


15.075 


15.054 


0.021 


0.1 


67.50 


15.071 


15.044 


0.027 


0.2 


40 Hz 


Test point 
locations [mm] 


^experimental \jYl/s] 


^theoretical l m / s ] 


Al/ [m/s] 


/^experimental % 


-67.50 


24.635 


24.647 


-0.012 


-0.0 


-54.00 


24.679 


24.655 


0.024 


0.1 


-22.50 


24.654 


24.663 


-0.009 


-0.0 


0.00 


24.667 


24.666 


0.001 


0.0 


22.50 


24.673 


24.663 


0.010 


0.0 


54.00 


24.686 


24.655 


0.031 


0.1 


67.50 


24.303 


24.647 


-0.344 


-1.4 


Table 6 Comparison between experimental and theoretical mean velocity values 




Vmean,th Ms] 


V 

v mean.exp 

[m/s] 


\AV\ 


\AV%\ 


10 Hz 


5.07 


4.95 


0.12 


2.42 


25Hz 


15.05 


14.73 


0.32 


2.17 


40 Hz 


24.65 


24.24 


0.41 


1.69 



The error estimated appears to be less than 3%, 
so validating the hypothesis of obtaining the flow 
velocity for the estimation of the flow rate 
directly from a unique measurement at a chosen 
position even in the squared cross -section pipe as 
well as in circular cross-section ones. 

The estimated position has resulted equal to 
about 0.78 of the half transversal dimension of 
the pipe. It is presumable that similar empirical 
model could be set even for rectangular cross- 



section pipes provided that its aspect ratio would 
not differ so much from unity. 

5. Conclusions 

An innovative empirical model has been set up 
for evaluating the suitable position where to 
measure the velocity of turbulent flow allowing to 
directly calculate the volumetric flow rate flowing 
in pipes with squared cross-section. The solving 
formula arises from the analogy to the circular 
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cross sections, for which the standard EN 12599 

3 

estimate such position at -R, R being the radius 
of the section. 

The model has been validated by measurements 
taken in a small subsonic wind tunnel whose test 
section has been properly designed. A laser 
Doppler anemometer has been utilized to measure 
the local velocities in the test section over a grid 
of points suitably chosen. 

The comparisons between calculations and 
measurements have shown discrepancies of order 
of 3% at most for the mean velocity and then for 
the volumetric flow rate, if the uncertainties on 
the section dimensions are neglected. Such 
discrepancies are much lesser than the errors 
accepted in common practice. 
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Abstract 

Strength measurements on the lower limb 
muscles are popular in the medical practice to 
evaluate the health status and effectiveness of 
training programs. Quality of strength measures 
obtained by Hand Held Dynamometer (HDD) is 
anyway under discussion, as reliability depends 
on operator as well as testing conditions. In this 
work the authors used an optoelectronic system to 
assist HHD measurements in order to check 
quality and to propose a method to reduce 
interval of uncertainties due to HHD method. 

1. Introduction 

Very often, in the medical practice and in 
rehabilitation programs, strength measurements 
are performed to evaluate the health status of 
patients and effectiveness of training programs. 
The force measurements by Hand Held 
Dynamometer (HHD) are also used to assess 
muscle strength in children with developmental 
disabilities, such as cerebral palsy (CP) [1]. 

The use of HHD is very popular because the 
device is inexpensive, easy to use and does not 
require specific patient preparation; in fact, the 
patient has just to assume the position defined by 
the protocol, e.g. sitting on a bench for knee 
flexion/extension trials, while the therapist applies 
the HHD on the leg and asks the patient to push 
against the HHD. At the end the therapist simply 
reads the force value on the display of the device. 
On the other side, HHD measures are affected by 
errors due to the operator and patient positioning. 

Studies were already conducted on the inter- 
tester reliability of the method and concluded that 
the method is questionable since due to the low 
reproducibility among trial repetitions [1], [2]. 

The reproducibility increases when the test is 
conducted by a trained clinician, but quality of 
measurements depends on application and 
positioning of the HHD [2] . 

The aim of this work is to conduct strength 
measurements while tracking position of the HHD 
with respect to the limb. Position and 3D 



orientation are recorded by an optoelectronic 
system. This will allow to quantify the error due 
to positioning, orientation and rotation of the 
HHD and to develop a method to correct the 
measurement. 

2. Methods 

2.1. Equipment 

Measures took place in the Motion Analysis and 
Robotics Laboratory at Ospedale Pediatrico 
Bambino Gesu, where a Vicon MX (Oxford, UK) 
motion capture system (8-camera-workstation, 
Nexus 1.7 software, 200 Hz, PluglnGait marker 
set based on the Davis protocol) is installed. 

A MicroFet™ dynamometer (Hoggan Scientific, 
Salt Lake City, UT) was equipped with four 
passive markers as shown in Fig. 1. Markers were 
placed on sticks fixed to the HHD in order to 
improve the visibility by the Vicon system. 

The central marker was used only for the static 
trial, to find the sensing HDD axis, and it was 
removed for the dynamic trials. 




Fig. 1: HHD Dynamometer equipped with passive 
markers. 

The Vicon system allowed to reconstruct the 
position of the subject and the HDD (Fig. 2) with 
an overall inaccuracy of 0.5 mm. Thus, the 
method proposed here allowed to compute the 
direction, orientation and application axis of the 
force measured by the HDD. 
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4. Conclusion 



Two adult healthy subjects were recruited for the 
tests. Knee extension trials were recorded with the 
subject sitting on a medical bench. Knee 
extension force was measured by applying the 
HHD on the shank, at ~ 50 mm from the ankle 
(Fig. 2). 




Fig. 2: Measured force, lever arm and knee torque. 

The subjects were requested to extend the knee 
against the dynamometer with as much force as 
they could, for ~5 s. The therapist had to push 
back in order to keep the shank still [1]. The test 
was repeated three times. Maximum force 
measured was recorded. The nominal moment 
was computed by multiplying the maximum force 
and the length of the shank, obtaining a constant 
value. Instead, the actual applied moment was 
computed by reconstructing the orientation of the 
HDD with respect to the subject. Knee extension 
moment was computed along the three anatomical 
axes: flex/extension, ab/adduction and intra/extra 
rotation obtaining a track for each axis, 
representing how the component changed during 
the trial. 

The angular range of motion of the knee was 
also measured, to ensure that the knee maintained 
the requested position during the trial. 

3. Results 

The time -histories of force and torques, 
determined in one of the performed trial, are 
depicted in Fig. 3 and 4. As expected, the main 
component of the force was in the antero- 
posterior direction (x-axis). The other components 
were lower (-1.5% and -30%). For the knee 
moment, the main component was along the z- 
axis (flex/ext -40 Nm). Other components were 
lower (-0% and -30%). The dotted line in Fig. 4 
is the nominal moment (-43 Nm). The black line 
is the intensity of moment vector. 



The preliminary results showed that the average 
torque in the principal direction was compatible 
with the respective nominal value. The difference 
was -3 Nm. The other components were lower. In 
case of erroneous application of HHD, the 
undesired components would increase. 

Further study and trial acquisition is required to 
define threshold values that may represent a 
numerical index to characterize the quality of a 
strength measurement. 

TRIAL 1 - Measured force expressed in CS rsh 
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Fig. 3: Subject No.l, test: right knee ext, repetition 
1: Components of measured force vector. 

TRIAL 1 - Knee moment expressed in CS rgh 







m knee,y 

Iml 


. . . — 





0 0.5 1 1.5 2 2.5 



time (s) 

Fig. 4: Subject No.l, test: right knee ext, repetition 
1: components of knee moment vector. 
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Abstract 

A Multicomponent Force Transducer (MFT) was 
devised by INRiM for a specific request in railway 
industry of measuring two force components 
(transversal and axial force). The measurement of 
the two transversal components was not required, 
but it was helpful for the practical use. INRiM, 
basing on its experience, developed and 
calibrated a full six-components hexapod-design 
prototype of MFT. Linear regression analysis was 
applied to get an empirical mathematical model 
linking the applied components to the MFT 
outputs. Best Subset Regression was adopted to 
identify parsimonious models by considering all 
possible combinations of predictor variables. 
Finally, the uncertainty associated to the values of 
transversal and axial force was evaluated. 

1. Introduction 

The present work was originated by a specific 
industrial request, i.e. the determination of 
"Chasse" (orthogonality error under load) 
generated by the transversal force due to 
compression on springs used in railway carriage 
suspension. The target uncertainties required by 
this application are 1% on the measurement of 
transversal force and class 1 of ISO 7500-1:2004 
[1] for axial force. Typical capacity being up to a 
few hundreds of kilonewtons, a Multicomponent 
Force Transducer (MFT) was devised ad hoc, on 
the basis of sizable experience in this field 
acquired at INRiM by developing several 
multicomponent measurement and calibration 
systems [2,3]. MFTs developed at INRiM found 
application in the characterization of primary 
force standard machines, e.g. control of parasitic 
components [4], as well as in such areas as 
robotics [5] and verification of material testing 
machines [6]. Both range (from a few newtons to 
hundreds of kilonewtons) and type (compound 
elements or monolithic) of MFTs developed cover 
a broad spectrum of applications. 

While measurement of main axial and 
transversal force magnitudes were mainly 
required, the direction of transversal force had 



also to be estimated, the direction of main axial 
force being the longitudinal axis of the spring. A 
three-component force transducer was therefore 
required. In addition, advantages coming from an 
even rough information on the applied moments 
allow a better identification of the calibration 
equations. Therefore, a full six-components 
hexapod-design prototype of MFT was 
developed. This design is typical of six degrees of 
freedom displacement measuring devices (Stewart 
platforms). It has been used till now in the field of 
MFTs mainly for low-force transducers such as 
those for application in robotics [7,8,9]; for the 
specific application an high load capacity MFT 
was developed in cooperation between INRiM 
and a company in the field of mechanical testing. 

2. Description of the MFT 

The MFT has an axial force capacity of 200 kN 
and of 30 kN for transversal force. Its built-up 
design is made up by six Uniaxial Force 
Transducers (UFTs) decoupled by elastic hinges 
(Fig. 1), substantially eliminating spurious 
components, which might otherwise affect single 
UFT measurements [10]. 




Fig. 1: Two of the six Uniaxial Force Transducers 
(UFTs) decoupled by elastic hinges at both ends. 



Such a structure (Fig. 2) enables measuring three 
force components (transversal, F x and F y , and 
axial, F z ), and three moment components, 
(bending, M x and M y , and torsion, M z ). 
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Fig. 2: Layout of the hexapod-shaped 
Multicomponent Force Transducer (MFT). 



A two-step spring test procedure applies first 
axial force, F z , and measures the two components 
of the transversal force (F x and F y ), enabling 
estimation of transversal force direction. The 
second step implies rotation of the spring in order 
to align transversal force with the x-axis, and 
applications of the axial force, iterating as 
required the rotation until a correct alignment is 
obtained (typically, no more than two iterations). 

3. Calibration of the MFT 

An important aspect which was considered for 
the choice of type of transducer design was the 
possibility of calibration, both initial and for 
maintenance over time. 

The main advantage of this type of MFT is that 
all the UFTs are subject to a tensile stress, greater 
or lower depending on the applied forces and 
moments. The first calibration is useful to 
determine both the effects of the geometry of the 
structure and the sensitivities of the UFTs. Since 
the MFT, installed in a spring testing machine 
(Fig. 3), is not moved over time, the geometry is 
not subject to significant changes. Thus, the 
subsequent calibrations can be restricted only to 
the evaluation of the variation over time of 
sensitivities of the UFTs by applying known 
values of vertical components F z . 

INRiM was tasked with the first calibration of 
the MFT in order to have traceable measurements, 
i.e. with an uncertainty evaluation [11]. A well- 
established metrological approach [2] asserts that 
known forces must be applied both independently 
and in combination to assess cross sensitivity 
among output channels (if any). Testing 
procedure strictly requires only calibration for F x 
and F z . 

Calibration for F z was performed on the primary 
INRiM deadweights standard machine with a 



capacity of 1 MN, following the international 
calibration procedure; results showed the MFT 
being in class 00 according to ISO 376:201 1 [12]. 
Calibration of F x was instead performed right on 
the spring testing machine, evaluating the all- 
important cross sensitivity with F z . The latter was 
generated using deadweights, transversal force 
being applied by mechanical devices and 
measured with a calibrated UFT (Fig. 3). 




Fig. 3: Calibration set-up on the spring testing 



machine. 



4. Analysis of calibration results 

Using the described calibration set-up (Fig. 3), 
transversal forces and axial forces were applied. 
Given the target of uncertainties (mentioned in 
section 1), and in the light of results obtained in 
preliminary tests, an abridged experimental plan 
with only two levels of F z and few levels of F x 
was deemed adequate and performed accordingly. 

According to the availability of deadweights in 
the site of calibration, the selected levels of F z 
were about 14 kN and 23 kN. With these levels of 
force, values of F x ranging up to 3 kN and 5 kN 
respectively have been chosen in order to apply 
only tensional stress to the single UFTs. By 
applying component F x an associated bending 
moment M y is also generated. The zero condition 
is obtained with the MFT in the minimum 
preloading condition. Applied components, 
corresponding force outputs and moment outputs 
are shown, respectively, in Table 1, 2 and 3. 
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Table 1. Applied components. 



n. 


.Fjc/kN 


F z /kN 


H AT /XT 

M v /J\rm 


i 
1 


A AAA 
U.UUU 


A AAA 
U.UUU 


A A 

U.U 


z 


A AAA 
U.UUU 


-zz.9o9 


A A 
U.U 


3 


A ^£.1 
O.JOJ 


OO A£A 


-zol.o 


A 

4 


1 ACA 

I.UjU 


TO A/CA 


CO /I o 
-JZ4.6 


J 


o a 1 a 


OO A£A 


1 AAA C 

-luuy.j 


0 


4.U4J 


00 AAA 


-ZUZZ.4 


7 


4.992 


-22.969 


-2495.9 


8 


0.000 


-14.066 


0.0 


9 


0.528 


-14.066 


-264.1 


10 


0.996 


-14.066 


-498.2 


11 


2.028 


-14.066 


-1014.2 


12 


3.006 


-14.066 


-1503.0 



Table 2. Force outputs corresponding to the applied 
components shown in Table 1. 



n. 


£/ x /kN 


Uy/W 


C/ z /kN 


1 


0.000 


0.000 


0.000 


2 


0.000 


0.000 


-22.969 


3 


0.520 


0.001 


-22.959 


4 


0.962 


0.004 


-22.959 


5 


1.850 


0.016 


-22.952 


6 


3.699 


0.049 


-22.931 


7 


4.563 


0.070 


-22.937 


8 


0.000 


0.000 


-14.066 


9 


0.483 


0.041 


-14.134 


10 


0.914 


0.044 


-14.134 


11 


1.862 


0.059 


-14.121 


12 


2.759 


0.075 


-14.111 



Table 3. Moment outputs corresponding to the 
applied components shown in Table 1. 



n. 


tWN-m 


U my /N-m 


U mz /N-m 


1 


0.0 


0.0 


0.0 


2 


0.0 


0.0 


0.0 


3 


76.7 


-279.4 


-0.4 


4 


69.7 


-517.7 


0.4 


5 


52.9 


-995.7 


-0.4 


6 


4.5 


-1996.3 


-0.4 


7 


-33.0 


-2465.2 


-4.7 


8 


0.0 


0.0 


0.0 


9 


-11.7 


-244.6 


7.7 


10 


-16.3 


-474.5 


8.4 


11 


-38.7 


-982.4 


7.7 


12 


-62.8 


-1462.2 


5.1 



While F x is expected to be related mainly to the 
output U x , and F z to the output U z , owing to cross- 
sensitivity the presence of other terms was 
anticipated. Linear regression analysis [13] was 
applied to get an empirical mathematical model 
linking the applied components to the MFT 
outputs. In particular, Best Subset Regression [14] 



was adopted to identify parsimonious models for 
F x and F z by considering all possible 
combinations of predictor variables (U Xi U y , U z , 
U mx , U my and U mz ). The limited number of 
experiments allowed only to consider models up 
to first order with interactions. This is acceptable 
considering the linearity characteristic of the 
UFTs. For the transversal force, it is obtained: 

F x =a-U x +b-U mx +c-U x -U z Eq. 1 

where a = 1.08, b = -6.16xl0" 2 m" 1 and 
c = -6.04xlO- 4 kN _1 . 
Instead, for the axial force, it is obtained: 

F z =d-U x +e-U y +f-U z Eq.2 

where d = -3.90x1 0 -2 , e = 2.00 and/= 1.00. 

The models include as significant terms for F x 
besides U x , also U mx and the interaction term 
U X U Z (Eq. 1), while the corresponding model for 
F z includes, besides U z , also U x and U y (Eq. 2). 
The significant contribution of the interaction 
U X U Z in Eq. 1 underlines the advantages of 
exploiting a full six-components MFT even when 
measurement of only two components is required. 

Once estimated the mathematical models, it was 
evaluated the expanded uncertainty associated to 
the values of F x and F z [15,16]. For all the 
equation coefficients, the relevant standard 
uncertainty was evaluated as the corresponding 
standard deviation obtained in the performed 
linear regression (Table 4). For the MFT outputs, 
the relevant standard uncertainty was evaluated 
considering the metrological characteristics of the 
UFTs and the geometry of the MFT (Table 5). 



Table 4. Standard deviations of equation 
coefficients. 



Coefficient 


Std. dev. 


a 


7.56xl0" 4 


b 


6.73xl0- 3 in 1 


c 


3.52xlO- 5 kN _1 


d 


2.65xl0 -3 


e 


1.19X10 -1 


f 


1.82xl0" 4 



Table 5. Relative standard uncertainties of MFT 
outputs. 



MFT 


Relative 


output 


std. unc. 


u x 


0.5% 




0.5% 


u z 


0.025% 


U mx 


5.9% 




2.9% 


u mz 


2.1% 
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The uncertainty evaluation according to 
GUM [11] and PUMA method (described in 
ISO 14253-2:2011 [17]) may be properly 
organized in a tabular format, with reference to 
EA-4/02 M:2013 [18]. A small modification from 
this format has been introduced by substituting 
standard deviations with variances; it has the 
advantage to manage additive quantities which 
can be compared more easily. In this way, 
individual contributions to variance of output 
quantity F x are shown for a specific working 
condition in Table 6. 



Symbols of independent variables appearing in 
the mathematical model and their values are 
written down in column Xj. Entries in column u(xj) 
are the standard uncertainties for each 
contribution, while values in column Vj represent 
the degrees of freedom. Coefficients of sensitivity 
Cj may be evaluated either by partial derivation, or 
numerically, and eventually contributions u/(F x ) 
of variance of dependent variable F x can be 
calculated. By taking into account all these 
information, it is possible to get the expanded 
uncertainty U(F X ). 



Table 6. Uncertainty table for the transversal force F x (expressed in kilonewton). 





u(xj) 


C J 


uf(F x ) 




u/(F x )/ vj 


Symbol 


Value 


a 


1.08 


7.6xl0" 4 


3.0 


5.1xl0" 6 


9 


2.9xl0" 12 


b 


-6.16xl0" 2 


6.7xl0" 3 


l.OxlO" 1 


4.5xl0" 7 


9 


2.3xl0" 14 


c 


-6.04xl0" 4 


3.5xl0" 5 


-42.0 


2.2xl0" 6 


9 


5.3xl0" 13 


u x 


3.00 


1.5xl0" 2 


1.1 


2.7xl0" 4 


61 


1.2xl0" 9 


u mx 


l.OOxlO" 1 


5.9xl0" 3 


-6.2xl0" 2 


1.3xl0" 7 


40 


4.5xl0" 16 


u z 


-14.00 


3.5xl0" 3 


-1.8xl0" 3 


4.1xl0" n 


100 


1.7xl0" 23 


F x 


3.26 




u 2 (F x ) 


2.7xl0" 4 




1.2xl0" 9 








u(F x ) 


1.7xl0" 2 


VFx 


64 


P 


95% 




tp( VFx) 


2.0 


U(F X ) 


3.3xl0" 2 


W(F X ) 


1.0% 



The same type of calculation was performed for 
F z . Therefore, it is obtained a relative expanded 
uncertainty of 1.0% for F x and 0.15% for F z . 
These values resulted to be nearly constant in all 
the working conditions. 

5. Conclusion 

A Multicomponent Force Transducer (MFT) was 
devised by INRiM for a specific request in 
railway industry. Even if the testing procedure 
strictly required only calibration for F x and F z , the 
INRiM experience suggested to develop a full six- 
components hexapod-design prototype of MFT. 
In fact, the analysis of calibration results showed 
that the measurement of F x is not only affected by 
the value of U x , but also by U mx and the 
interaction term U X U Z . 

Besides, the uncertainty associated to the values 
of transversal and axial force was evaluated. The 
applied calibration method allowed to keep these 
uncertainty values within the targets required by 
the customer. In particular, it is obtained a relative 
expanded uncertainty of 1.0% for transversal 



force and 0.15% for axial force, much better than 
the required class 1 of ISO 7500-1:2004. 
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